- Email: [email protected]
Free energy profiles of cocaine esterase-cocaine binding process by molecular dynamics and potential of mean force simulations
Chemico-Biological Interactions 259 (2016) 142e147
Contents lists available at ScienceDirect
Chemico-Biological Interactions journal homepage: www.elsevier.com/locate/chembioint
Free energy profiles of cocaine esterase-cocaine binding process by molecular dynamics and potential of mean force simulations Yuxin Zhang a, b, c, 1, Xiaoqin Huang b, 1, Keli Han a, Fang Zheng b, c, Chang-Guo Zhan b, c, * a
State Key Laboratory of Molecular Reaction Dynamics, Dalian Institute of Chemical Physics, Chinese Academy of Science, Zhongshan Road 457, Dalian 116023, PR China b Molecular Modeling and Biopharmaceutical Center, College of Pharmacy, University of Kentucky, 789 South Limestone Street, Lexington, KY 40536, USA c Department of Pharmaceutical Sciences, College of Pharmacy, University of Kentucky, 789 South Limestone Street, Lexington, KY 40536, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 9 January 2016 Received in revised form 29 April 2016 Accepted 5 May 2016 Available online 6 May 2016
The combined molecular dynamics (MD) and potential of mean force (PMF) simulations have been performed to determine the free energy profile of the CocE)-(þ)-cocaine binding process in comparison with that of the corresponding CocE-()-cocaine binding process. According to the MD simulations, the equilibrium CocE-(þ)-cocaine binding mode is similar to the CocE-()-cocaine binding mode. However, based on the simulated free energy profiles, a significant free energy barrier (~5 kcal/mol) exists in the CocE-(þ)-cocaine binding process whereas no obvious free energy barrier exists in the CocE-()-cocaine binding process, although the free energy barrier of ~5 kcal/mol is not high enough to really slow down the CocE-(þ)-cocaine binding process. In addition, the obtained free energy profiles also demonstrate that (þ)-cocaine and ()-cocaine have very close binding free energies with CocE, with a negligible difference (~0.2 kcal/mol), which is qualitatively consistent with the nearly same experimental KM values of the CocE enzyme for (þ)-cocaine and ()-cocaine. The consistency between the computational results and available experimental data suggests that the mechanistic insights obtained from this study are reasonable. © 2016 Elsevier Ireland Ltd. All rights reserved.
Keywords: Esterase Cocaine Binding process Energy barrier
1. Introduction Cocaine may have two enantiomers, ()-cocaine and (þ)-cocaine (see Fig. 1). The naturally occurring and biologically active ()-cocaine is known as the most addictive substance abused by millions of people worldwide [1e4]. Despite of decades of research efforts, there is no FDA-approved anti-cocaine therapeutic agent specific for cocaine addiction or overdose [4,5]. As well known, cocaine exerts its effects on central nervous system (CNS) through blocking the reuptake of neurotransmitter dopamine, thus potentiating the effects of dopamine in the synapse [6]. Traditional pharmacodynamic approaches have failed to produce a therapeutically useful small-molecule drug because of the difficulties inherent in blocking a blocker like cocaine. As an alternative,
* Corresponding author. Department of Pharmaceutical Sciences, College of Pharmacy, University of Kentucky, 789 South Limestone Street, Lexington, KY 40536, USA. E-mail address: [email protected] (C.-G. Zhan). 1 These authors contributed equally to this work. http://dx.doi.org/10.1016/j.cbi.2016.05.011 0009-2797/© 2016 Elsevier Ireland Ltd. All rights reserved.
pharmacokinetic approach using an efficient cocaine-metabolizing enzyme has been recognized as a promising strategy for therapeutic treatment of cocaine overdose and abuse. The enzyme strategy aims at accelerating cocaine metabolism (particularly through catalyzing cocaine benzoyl ester hydrolysis) and, therefore, eliminating cocaine rapidly in the plasma [1,7e9]. For this purpose, cocaine esterase (CocE) [10e13] is a promising choice as a potential anti-cocaine agent for therapeutic treatment of cocaine overdose and abuse, because CocE is the most efficient natural enzyme against ()-cocaine [12]. Native CocE is a bacterial enzyme with a catalytic center similar to that of butyrylcholinesterase (BChE) [14,15], but the catalytic efficiency of CocE [12] is two orders of magnitude higher than that of wild-type BChE against ()-cocaine [16]. Native CocE and the designed thermostable mutants are capable of protecting against cocaine-induced lethality [12,13,17,18]. Thus, both BChE and CocE are interesting protein scaffolds that may be used to design enzyme mutants with improved catalytic efficiency against ()-cocaine. To rationally redesign an enzyme (BChE or CocE) for development of anti-cocaine therapeutics, one first needs to understand the detailed catalytic mechanisms of the enzyme against both
Y. Zhang et al. / Chemico-Biological Interactions 259 (2016) 142e147
143
2. Methods 2.1. MD simulation
Fig. 1. Molecular structures of cocaine enantiomers.
()-cocaine and (þ)-cocaine and understand their similarities and differences in the enzyme-substrate binding and subsequent transformation processes. Previously reported studies [7,14,15,19,20] revealed that BChE- and CocE-catalyzed hydrolyses of ()-cocaine and (þ)-cocaine share a similar catalytic reaction pathway for the chemical reaction process. For each of the enzymatic reactions, the chemical reaction process always includes both the acylation and deacylation stages, regardless of the enzyme (CocE or BChE) and substrate (()-cocaine or (þ)-cocaine). There is also no significant difference in Michaelis-Menten constant (KM) between BChE-catalyzed ()-cocaine hydrolysis and BChEcatalyzed (þ)-cocaine hydrolysis [21]. However, the enzymesubstrate binding processes are remarkably different [22,23]. For example, it has been known that the enzyme-substrate binding process is the rate-determining step for BChE-catalyzed hydrolysis of ()-cocaine, whereas the chemical reaction process is the ratedetermining step for BChE-catalyzed hydrolysis of (þ)-cocaine [16,24,25]. Based on the understanding of the mechanistic difference, we were able to perform structure and mechanism based BChE redesign, and successfully designed and discovered new BChE mutants with considerably improved catalytic efficiency against ()-cocaine [25e32]. It has also been known that the CocE-()-cocaine binding process [33] is remarkably different from the BChE-()-cocaine binding process [23]. BChE-()-cocaine binding process is associated with a very high energy barrier [22,23], whereas the CocE()-cocaine binding process is barrier-less [33]. For the remaining mechanistic question to be addressed, it is unknown whether the CocE-(þ)-cocaine binding process is also different from the CocE()-cocaine binding process, because no computational study has been reported on the CocE-(þ)-cocaine binding process. Understanding the difference/similarity between the CocE-(þ)-cocaine and CocE-()-cocaine binding processes would be valuable for guiding further anti-cocaine enzyme redesign efforts. In the present study, we carried out molecular dynamics (MD) and potential of mean force (PMF) simulations to determine the free energy profile of the CocE-(þ)-cocaine binding process in comparison with the corresponding CocE-()-cocaine binding free energy profile. The combined MD and PMF simulations have revealed some interesting difference and similarity in the free energy profile between the CocE-(þ)-cocaine and CocE-()-cocaine binding processes.
The CocE structure used in this study is a thermally stable CocE mutant (T172R/G173Q) designed in our previous study [12]. It has been known that the T172R/G173Q mutations can significantly stabilize the CocE structure without changing the catalytic activity. The same T172R/G173Q mutant CocE was also used to study the CocE-()-cocaine binding process [33]. For convenience, below we simply refer CocE to the T172R/G173Q mutant CocE, unless indicated explicitly otherwise. Initial structure of the CocE-(þ)-cocaine binding complex was prepared based on the available X-ray crystal structure [34] of CocE (PDB code 3I2F for the T172R/G173Q mutant CocE) and the results of our previous molecular docking and MD simulations on CocE binding with ()-cocaine [14,33]. The initial CocE-(þ)-cocaine binding structure was energy-minimized by using the Sander module of Amber 12 program [35] via a combined use of the steepest descent/conjugate gradient algorithms, with a convergence criterion of 0.01 kcal mol1 Å1, and the non-bonded cutoff distance was set to 10.0 Å. MD simulation was performed by using the Sander module of the Amber 12 program package. The general procedure of the MD simulation was similar to that used in our previously reported other computational studies [12,13,33,36]. In particular, the molecular mechanics force field parameters and the partial charges of (þ)-cocaine atoms were adopted directly from those developed in our previous studies on other complexes [12,13,33,36]. Briefly, the partial charges of (þ)-cocaine atoms were calculated by using the restrained electrostatic potential-fitting (RESP) protocol implemented in the Antechamber module of the Amber 12 program, following the electrostatic potential (ESP) calculation at ab initio HF/6-31G* level using Gaussian 03 program [37]. For the energy minimization and MD simulation, the CocE(þ)-cocaine complex was solvated in an orthorhombic box of TIP3P water molecules [38] with a minimum solvent-wall distance of 10 Å. Sodium counter ions (Naþ) were added to neutralize the solvated system. The solvated system was gradually heated to 298.15 K by the weak-coupling method [39] and equilibrated for 400 ps. During the MD simulations, a 10.0 Å non-bonded interaction cutoff was used and the non-bonded list was updated every 25 steps. The motion for the mass center of the system was removed every 1000 steps. The particle-mesh Ewald (PME) method [40,41] was used to treat long-range electrostatic interactions. The lengths of covalent bonds involving hydrogen atoms were fixed with the SHAKE algorithm [42], enabling the use of a 2-fs time step to numerically integrate the equations of motion. Finally, the production MD was kept running for ~2.0 ns with a periodic boundary condition in the NTP ensemble at T ¼ 298.15 K with Berendsen temperature coupling, and at P ¼ 1 atm with isotropic moleculebased scaling [40,43]. 2.2. Potential of mean force (PMF) simulation In order to explore the free energy profiles for the process of (þ)-cocaine binding with CocE, the PMF simulation was carried out by using umbrella-sampling [44] MD simulation. The classic PMF definition [45] can be represented by a function of reaction coordinate as
uðcÞ ¼ RT ln〈rðcÞ〉 UðcÞ þ F
(1)
in which r(c) is the probability density along the reaction coordinate c, R is the gas constant, T is the simulation temperature, U(c) is the biasing potential applied in the umbrella-sampling MD simulation, and F is the normalization constant. According to this
144
Y. Zhang et al. / Chemico-Biological Interactions 259 (2016) 142e147
approach, the reaction coordinate is usually divided into different regions, i.e., windows, and each window is sampled separately. A biasing (umbrella) potential, i.e. U(c), is applied for each window in order to obtain nearly uniform sampling of the potential energy surface. In the present study, the reaction coordinate was defined as the distance from the mass center of the non-hydrogen atoms of (þ)-cocaine to the mass center of the non-hydrogen atoms on the side chains of residues H87, V121, and L146 of CocE. The total number of windows was 73, with a window size of 0.3 Å, covering the reaction coordinate values from ~13.385 Å to ~34.985 Å. The biasing force constant applied in different windows of umbrellasampling was 10.0 kcal/(molÅ2), as used in our previous PMF simulations on the CocE-()-cocaine binding process [33]. For each umbrella-sampling window, the initial complex structure was selected from the last snapshot of the PMF simulations of the previous window. The selected structure for each window was first equilibrated for 200 ps and then kept running for 800 ps for production sampling. The frequency for data collection was set to 1 fs, which was the same as that of the time step of the umbrellasampling MD. Previous computational study [33] on the CocE()-cocaine binding process demonstrated that this PMF protocol led to converged PMF free energy profiles, as further longer MD and PMF simulations did not significantly change the obtained free energy profiles. After the umbrella-sampling MD simulations for all windows were completed, the data collected from separate simulation windows were combined along the reaction coordinate. These data were then used to calculate the PMF for the whole binding process with the weighed histogram analysis method (WHAM) [46,47] using the code developed by Alan Grossfield (http://membrane.urmc.rochester.edu/Software/WHAM/WHAM. html). Notably, the equation for the biasing potential implemented in the Amber program for the PMF simulations is U(c) ¼ kAmber (cc0)2, whereas the corresponding equation implemented in the WHAM program for the PMF free energy analysis is U(c) ¼ (kWHAM/2) (cc0)2. So, the equations implemented in the two programs are inconsistent in the force constant. Clearly, kWHAM ¼ 2kAmber. In this study, we used kAmber ¼ 10 kcal/ mol and kWHAM ¼ 20 kcal/mol. It should be noted that such a difference between kWHAM and kAmber was not noted in the previous analysis [33] of the PMF data obtained for the CocE-()-cocaine binding. As a result, the previous analysis [33] of the PMF data obtained for the CocE-()-cocaine binding was performed inconsistently by using kAmber ¼ 10 kcal/mol and kWHAM ¼ 10 kcal/mol. In this way, the previous analysis using the WHAM program systematically under-estimated the PMF free energy profile (and the binding free energy) for the CocE-()-cocaine binding by a factor of 2, although the systematical under-estimation did not change any of the conclusions in the previous study [33]. All of the MD and umbrella-sampling MD simulations were performed on a supercomputer (a DELL Cluster with 388 nodes or 4816 processors) at the University of Kentucky’s Computer Center.
Fig. 2. Tracked internuclear distances in the MD-simulated CocE-(þ)-cocaine binding complex. (A) C32(COC)/S117 (black color) represents the distance between the carbonyl carbon on the benzoyl group of (þ)-cocaine and the hydroxyl oxygen of the S117 side chain; O33(COC)/Y44 (red color) refers to the distance between the carbonyl oxygen on the benzoyl group of (þ)-cocaine and the hydroxyl hydrogen of the Y44 side chain; O33(COC)/Y118 (blue color) represents the distance between the carbonyl oxygen on the benzoyl group of (þ)-cocaine and the backbone hydrogen of Y118. (B) S117/H287 (black color) and H287—D259 (red color) represent the distances related to hydrogen bonds within the catalytic triad S117-H287-D259. Y44/W166 (blue color) represents the hydrogen bond distance between the hydroxyl oxygen on the Y44 side chain and the side chain of W166. The block colored distances are not visible because they overlap with the blue colored distances. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3. Results and discussion 3.1. CocE-(þ)-cocaine binding structure According to the MD-simulated CocE-(þ)-cocaine binding structure, (þ)-cocaine binding site is located on the interface of the three domains of CocE (domains I to III). Depicted in Fig. 2 are important distances in the simulated CocE-(þ)-cocaine complex, and depicted in Fig. 3 is the structure in the last snapshot of the MD simulated CocE-(þ)-cocaine complex, revealing the protein environment of (þ)-cocaine in the complex. Figs. 2 and 3 clearly show the important hydrogen-bonding interactions that were very stable
Fig. 3. MD-simulated CocE-(þ)-cocaine binding structure. The enzyme is shown as golden ribbons, and (þ)-cocaine is in ball-and-stick style. Key residues of the enzyme within 5 Å around the (þ)-cocaine molecule are shown in stick style and colored by atom types. Hydrogen bonds are represented as dashed lines with labeled distances (in unit of Å). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Y. Zhang et al. / Chemico-Biological Interactions 259 (2016) 142e147
during the MD simulation process. In the simulated CocE(þ)-cocaine complex (Fig. 3), the benzoyl ester group of (þ)-cocaine is located in a sub-binding site composed of hydrophobic residues W151, W166, L169, L407, F408, F261, and P150 of CocE, packing in parallel with the aromatic side chain of W166, and in perpendicular with the aromatic side chain of F261. As shown in Fig. 3, S117, H287, and D259 residues form a stable catalytic triad with a stable hydrogen bonding network, and the carbonyl carbon of (þ)-cocaine benzoyl ester group is about ~3.5 Å away from the hydroxyl oxygen of S117 side chain. ~3.5 Å is an ideal distance for the hydroxyl oxygen to initiate the nucleophilic attack on the carbonyl carbon, according to our previous computational studies on other ester hydrolysis reactions catalyzed by a serine esterase (with a similar catalytic triad) like BChE [7,24]. In addition, the carbonyl oxygen on the benzoyl ester group of (þ)-cocaine interacts with the hydroxyl group of Y44 side chain through hydrogen-bonding interactions. As seen in Figs. 2 and 3, there is a strong hydrogen bond between the carbonyl oxygen of (þ)-cocaine benzoyl ester group and the hydroxyl group of Y44 side chain, with an O … H distance of ~1.9 Å. This hydrogen bond is expected to become stronger and stronger in going from the prereaction CocE-()-cocaine binding to the transition state and intermediate due to the increasing negative charge in the carbonyl oxygen atom during the acylation process. Notably, the CocE(þ)-cocaine binding mode simulated in this study is generally similar to the corresponding CocE-()-cocaine binding mode found in our previous study [33]. 3.2. Free energy profile of the CocE-(þ)-cocaine binding process In order to calculate the free energy of the CocE-(þ)-cocaine binding, the PMF simulations were carried out starting from the MD-simulated CocE-(þ)-cocaine complex structure. Based on the data collected from the umbrella-sampling MD simulations, the PMF for the CocE-(þ)-cocaine binding was determined. The distance between the mass center of the non-hydrogen atoms of (þ)-cocaine and the mass center of the non-hydrogen atoms on the side chains of residues H87, V121, and L146 of the enzyme was used as the reaction coordinate for the PMF calculations. The selection of the reaction coordinate was based on the structural features of our modeled CocE-(þ)-cocaine binding structures in the present study, because we found that the direction of reaction coordinate roughly went through the central point of the active site of CocE to reach the molecular surface of the enzyme. For a fair comparison between the free energy profile of the CocE-(þ)-cocaine binding process and that of the CocE-()-cocaine binding process, we also analyzed the previously obtained PMF data (based on the use of the same computational protocol described in the present study) for the CocE-()-cocaine binding process with kWHAM ¼ 20 kcal/mol (because the PMF simulations [33] were carried out using kAmber ¼ 10 kcal/mol). Depicted in Fig. 4 are the PMF free energy profiles obtained for the CocE-(þ)-cocaine and CocE-()-cocaine binding processes. Based on the PMF free energy profiles depicted in Fig. 4, there is a significant free energy barrier (~2 kcal/mol) when the reaction coordinate c ¼ ~17.6 Å (associated with a transition state) in the CocE-(þ)-cocaine binding process, whereas we cannot identify an obvious free energy barrier in the CocE-()-cocaine binding process. The free energy barrier for the CocE-(þ)-cocaine binding process is associated with (þ)-cocaine rotation in the binding pocket. Fig. 5(A) depicts the CocE-(þ)-cocaine complex structure when c ¼ ~17.8 Å (right before the transition state), and Fig. 5(B) shows the CocE-(þ)-cocaine complex structure when c ¼ ~17.3 Å (right after the transition state). Further, in terms of the binding kinetics, a free energy barrier of
145
Fig. 4. Simulated free energy profile of the CocE-(þ)-cocaine binding process (A) in comparison with that of the CocE-()-cocaine binding process (B). The reaction coordinate was defined as the distance between the mass center of the non-hydrogen atoms of (þ)/()-cocaine and the mass center of the non-hydrogen atoms on the side chains of residues H87, V121, and L146 of the enzyme.
~2 kcal/mol is associated with a rate constant as large as ~3 1011 s1 (faster than the diffusion) according to the classical rate constant theory [48], suggesting that the binding process is much faster than the chemical reaction process (with a much lower catalytic rate constant kcat [10,12]). So, the binding process associated with a free energy barrier of ~2 kcal/mol is still very fast and, hence, easy to reach the thermodynamic equilibration of the CocE(þ)-cocaine binding. So, the free energy barrier of ~2 kcal/mol is not high enough to really slow down the enzymatic reaction process (compared to the much slower chemical reaction process), and the CocE-substrate binding process is not rate-determining for both (þ)-cocaine and ()-cocaine. In this sense, we may also say that the CocE-(þ)-cocaine and CocE-()-cocaine binding processes are similar, in contrast to the remarkable difference between BChE(þ)-cocaine and BChE-()-cocaine binding processes. It has been known that the BChE-()-cocaine binding process is ratedetermining, whereas the BChE-(þ)-cocaine binding process is not [22e24]. In terms of the equilibrium CocE-substrate binding, the free energy profiles depicted in Fig. 4 indicate that the calculated CocE-substrate binding free energy is 12.2 kcal/mol for (þ)-cocaine and 12.4 kcal/mol for ()-cocaine. The difference (~0.2 kcal/mol) between (þ)-cocaine and ()-cocaine is negligible, which is qualitatively consistent with the nearly same experimental KM values of the CocE enzyme for the two substrates (KM ¼ ~15 mM for (þ)-cocaine and KM ¼ ~13 mM for ()-cocaine) [20]. 4. Conclusion The combined MD and PMF simulations have led to determination of the free energy profile of the CocE-(þ)-cocaine binding process in comparison with that of the CocE-()-cocaine binding
146
Y. Zhang et al. / Chemico-Biological Interactions 259 (2016) 142e147
DA041115) and R01 grants (R01 DA035552, R01 DA032910, R01 DA013930, and R01 DA025100). The authors also acknowledge the Computer Center at University of Kentucky for supercomputing time on a Dell Supercomputer Cluster consisting of 388 nodes or 4816 processors. Transparency document Transparency document related to this article can be found online at http://dx.doi.org/10.1016/j.cbi.2016.05.011. References
Fig. 5. Simulated CocE-(þ)-cocaine binding structures right before and after the rotational transition state (when c ¼ ~17.6 Å): (A) the CocE-(þ)-cocaine binding structure when c ¼ ~17.8 Å (right before the transition state); (B) the CocE-(þ)-cocaine binding structure c ¼ ~17.3 Å (right after the transition state).
process. According to the obtained free energy profiles, a significant free energy barrier (~5 kcal/mol) exists in the CocE-(þ)-cocaine binding process, whereas no obvious free energy barrier exists in the CocE-()-cocaine binding process. The free energy barrier for the CocE-(þ)-cocaine binding process is associated with (þ)-cocaine rotation in the binding pocket. On the other hand, the free energy barrier of ~2 kcal/mol is not high enough to really slow down the CocE-(þ)-cocaine binding process, and the CocEsubstrate binding process should not be rate-determining for both (þ)-cocaine and ()-cocaine. The obtained free energy profiles also reveal that (þ)-cocaine and ()-cocaine have very close binding free energies with CocE. The negligible binding free energy difference (~0.2 kcal/mol) between (þ)-cocaine and ()-cocaine is qualitatively consistent with the nearly same experimental KM values of the CocE enzyme for (þ)-cocaine and ()-cocaine. The consistency between the computational and experimental data suggests that the mechanistic insights obtained from the computational simulations in this study are reasonable. Competing financial interests statement The authors declare that there is no conflict of interest for this work. Acknowledgments This work was supported by the National Institutes of Health (NIH) through the NIDA Translational Avant-Garde Award (UH2
[1] D.W. Landry, K. Zhao, G.X. Yang, M. Glickman, T.M. Georgiadis, Antibodycatalyzed degradation of cocaine, Science 259 (1993) 1899e1901. [2] A.L. Milton, B.J. Everitt, Wiping drug memories, Science 336 (2012) 167e168. [3] K.D. Ersche, P.S. Jones, G.B. Williams, A.J. Turton, T.W. Robbins, E.T. Bullmore, Abnormal brain structure implicated in stimulant drug addiction, Science 335 (2012) 601e604. [4] F. Zheng, C.-G. Zhan, Are pharmacokinetic approaches feasible for treatment of cocaine addiction and overdose? Future Med. Chem. 4 (2012) 125e128. [5] D.A. Gorelick, Pharmacokinetic approaches to treatment of drug addiction, Exp. Rev. Clin. Pharmacol. 1 (2008) 277e290. [6] R. Chen, M.R. Tilley, H. Wei, F. Zhou, F.-M. Zhou, S. Ching, N. Quan, R.L. Stephens, E.R. Hill, T. Nottoli, D.D. Han, H.H. Gu, Abolished cocaine reward in mice with a cocaine-insensitive dopamine transporter, Proc. Natl. Acad. Sci. U. S. A. 103 (2006) 9333e9338. [7] C.-G. Zhan, F. Zheng, D.W. Landry, Fundamental reaction mechanism for cocaine hydrolysis in human butyrylcholinesterase, J. Am. Chem. Soc. 125 (2003) 2462e2474. [8] F. Zheng, C.-G. Zhan, Modeling of pharmacokinetics of cocaine in human reveals the feasibility for development of enzyme therapies for drugs of abuse, PLoS Comput. Biol. 8 (2012) e1002610. [9] X. Chen, L. Xue, S. Hou, Z. Jin, T. Zhang, F. Zheng, C.-G. Zhan, Long-acting cocaine hydrolase for addiction therapy, Proc. Natl. Acad. Sci. U. S. A. 113 (2016) 422e427. [10] N.A. Larsen, J.M. Turner, J. Stevens, S.J. Rosser, A. Basran, R.A. Lerner, N.C. Bruce, I.A. Wilson, Crystal structure of a bacterial cocaine esterase, Nat. Struct. Mol. Biol. 9 (2002) 17e21. [11] M.C. Ko, L.D. Bowen, D. Narasimhan, A.A. Berlin, N.W. Lukacs, R.K. Sunahara, Z.D. Cooper, J.H. Woods, Cocaine esterase: interactions with cocaine and immune responses in mice, J. Pharmacol. Exp. Ther. 320 (2007) 926e933. [12] D. Gao, D.L. Narasimhan, J. Macdonald, M.-C. Ko, D.W. Landry, J.H. Woods, R.K. Sunahara, C.-G. Zhan, Thermostable variants of cocaine esterase for longtime protection against cocaine toxicity, Mol. Pharmacol. 75 (2009) 318e323. [13] L. Fang, K.M. Chow, S. Hou, L. Xue, D.W. Rodgers, F. Zheng, C.-G. Zhan, Rational design, preparation, and characterization of a therapeutic enzyme mutant with improved stability and function for cocaine detoxification, ACS Chem. Biol. 9 (2014) 1764e1772. [14] J.J. Liu, A. Hamza, C.-G. Zhan, Fundamental reaction mechanism free energy profile for (-)-cocaine hydrolysis catalyzed by cocaine esterase, J. Am. Chem. Soc. 131 (2009) 11964e11975. [15] J. Liu, C.-G. Zhan, Reaction pathway and free energy profile for cocaine hydrolase-catalyzed hydrolysis of (e)-cocaine, J. Chem. Theory Comput. 8 (2012) 1426e1435. [16] H. Sun, J. El Yazal, O. Lockridge, L.M. Schopfer, S. Brimijoin, Y.-P. Pang, Predicted michaelis-menten complexes of cocaine-butyrylcholinesterase: engineering effective butyrylcholinesterase mutants for cocaine detoxication, J. Biol. Chem. 276 (2001) 9330e9336. [17] G.T. Collins, R.L. Brim, D. Narasimhan, M.-C. Ko, R.K. Sunahara, C.-G. Zhan, J.H. Woods, Cocaine esterase prevents cocaine-induced toxicity and the ongoing intravenous self-administration of cocaine in rats, J. Pharmacol. Exp. Ther. 331 (2009) 445e455. [18] R.L. Brim, M.R. Nance, D.W. Youngstrom, D. Narasimhan, C.-G. Zhan, J.J. Tesmer, R.K. Sunahara, J.H. Woods, A thermally stable form of bacterial cocaine esterase: a potential therapeutic agent for treatment of cocaine abuse, Mol. Pharmacol. 77 (2010) 593e600. [19] C.-G. Zhan, D. Gao, Catalytic mechanism and energy barriers for butyrylcholinesterase-catalyzed hydrolysis of cocaine, Biophys. J. 89 (2005) 3863e3872. [20] J.J. Liu, X.Y. Zhao, W.C. Yang, C.-G. Zhan, Reaction mechanism for cocaine esterase-catalyzed hydrolyses of (þ)- and (-)-cocaine: unexpected common rate-determining step, J. Phys. Chem. B 115 (2011) 5017e5025. [21] L. Xue, S. Hou, W. Yang, L. Fang, F. Zheng, C.-G. Zhan, Catalytic activities of a cocaine hydrolase engineered from human butyrylcholinesterase against (þ)and (-)-cocaine, Chem. Biol. Interact. 203 (2013) 57e62. [22] X. Huang, Y. Pan, F. Zheng, C.-G. Zhan, Reaction pathway and free energy profile for pre-chemical reaction step of human butyrylcholinesterasecatalyzed hydrolysis of (-)-cocaine by combined targeted molecular dynamics and potential of mean force simulations, J. Phys. Chem. B 114 (2010) 13545e13554.
Y. Zhang et al. / Chemico-Biological Interactions 259 (2016) 142e147 [23] X. Huang, F. Zheng, C.-G. Zhan, Human butyrylcholinesterase-cocaine binding pathway and free energy profiles by molecular dynamics and potential of mean force simulations, J. Phys. Chem. B 115 (2011) 11254e11260. [24] A. Hamza, H. Cho, H.-H. Tai, C.-G. Zhan, Molecular dynamics simulation of cocaine binding with human butyrylcholinesterase and its mutants, J. Phys. Chem. B 109 (2005) 4776e4782. [25] Y. Pan, D. Gao, W. Yang, H. Cho, G. Yang, H.-H. Tai, C.G. Zhan, Computational redesign of human butyrylcholinesterase for anticocaine medication, Proc. Natl. Acad. Sci. U. S. A. 102 (2005) 16656e16661. [26] F. Zheng, W.C. Yang, M.C. Ko, J.J. Liu, H. Cho, D.Q. Gao, M. Tong, H.H. Tai, J.H. Woods, C.-G. Zhan, Most efficient cocaine hydrolase designed by virtual screening of transition states, J. Am. Chem. Soc. 130 (2008) 12148e12155. [27] F. Zheng, W. Yang, L. Xue, S. Hou, J. Liu, C.-G. Zhan, Design of high-activity mutants of human butyrylcholinesterase against (-)-cocaine: structural and energetic factors affecting the catalytic efficiency, Biochemistry 49 (2010) 9113e9119. [28] L. Xue, M.-C. Ko, M. Tong, W. Yang, S. Hou, L. Fang, J. Liu, F. Zheng, J.H. Woods, H.-H. Tai, C.-G. Zhan, Design, preparation, and characterization of high-activity mutants of human butyrylcholinesterase specific for detoxification of cocaine, Mol. Pharmacol. 79 (2011) 290e297. [29] W. Yang, Y. Pan, L. Fang, D. Gao, F. Zheng, C.-G. Zhan, Free-energy perturbation simulation on transition states and high-activity mutants of human butyrylcholinesterase for (-)-cocaine hydrolysis, J. Phys. Chem. B 114 (2010) 10889e10896. [30] W. Yang, Y. Pan, F. Zheng, H. Cho, H.-H. Tai, C.-G. Zhan, Free-energy perturbation simulation on transition states and redesign of butyrylcholinesterase, Biophys. J. 96 (2009) 1931e1938. [31] F. Zheng, L. Xue, S. Hou, J. Liu, M. Zhan, W. Yang, C.-G. Zhan, A highly efficient cocaine-detoxifying enzyme obtained by computational design, Nat. Commun. 5 (2014) 3457, http://dx.doi.org/10.1388/ncomms4457. [32] D. Gao, H. Cho, W. Yang, Y. Pan, G. Yang, H.H. Tai, C.-G. Zhan, Computational design of a human butyrylcholinesterase mutant for accelerating cocaine hydrolysis based on the transition-state simulation, Angew. Chem. Int. Ed. 45 (2006) 653e657. [33] X. Huang, X. Zhao, F. Zheng, C.-G. Zhan, Cocaine esterase-cocaine binding process and the free energy profiles by molecular dynamics and potential of mean force simulations, J. Phys. Chem. B 116 (2012) 3361e3368. [34] D. Narasimhan, M.R. Nance, D. Gao, M.-C. Ko, J. Macdonald, P. Tamburi, D. Yoon, D.M. Landry, J.H. Woods, C.-G. Zhan, J.J.G. Tesmer, R.K. Sunahara, Structural analysis of thermostabilizing mutations of cocaine esterase, Protein. Eng. Des. Sel. 23 (2010) 537e547. [35] D.A. Case, T.A. Darden, T.E. Cheatham Iii, C.L. Simmerling, J. Wang, R.E. Duke, R. Luo, R.C. Walker, W. Zhang, K.M. Merz, B. Roberts, S. Hayik, A. Roitberg, ry, K.F. Wong, F. Paesani, J. Vanicek, G. Seabra, J. Swails, A.W. Goetz, I. Kolossva R.M. Wolf, J. Liu, X. Wu, S.R. Brozell, T. Steinbrecher, H. Gohlke, Q. Cai, X. Ye, J. Wang, M.J. Hsieh, G. Cui, D.R. Roe, D.H. Mathews, M.G. Seetin, R. SalomonFerrer, C. Sagui, V. Babin, T. Luchko, S. Gusarov, A. Kovalenko, P.A. Kollman, AMBER 12, University of California, San Francisco, 2012.
147
[36] J. Liu, A. Hamza, C.-G. Zhan, Fundamental reaction mechanism free energy profile for ()-cocaine hydrolysis catalyzed by cocaine esterase, J. Am. Chem. Soc. 131 (2009) 11964e11975. [37] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, J.A. Montgomery Jr., T. Vreven, K.N. Kudin, J.C. Burant, J.M. Millam, S.S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G.A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J.E. Knox, H.P. Hratchian, J.B. Cross, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, P.Y. Ayala, K. Morokuma, G.A. Voth, P. Salvador, J.J. Dannenberg, V.G. Zakrzewski, S. Dapprich, A.D. Daniels, M.C. Strain, O. Farkas, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J.V. Ortiz, Q. Cui, A.G. Baboul, S. Clifford, J. Cioslowski, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, M.W. Wong, C. Gonzalez, J.A. Pople, Gaussian 03, Rev A.1, Gaussian Inc., Pittsburgh, PA, 2003. [38] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys. 79 (1983) 926e935. [39] T. Morishita, Fluctuation formulas in molecular dynamics simulations with the weak coupling heat bath, J. Chem. Phys. 113 (2000) 2976e2982. [40] T. Darden, D. York, L. Pedersen, Particle mesh ewalddan N_log(N) method for ewald sums in large systems, J. Chem. Phys. 98 (1993) 10089e10092. [41] A. Toukmaji, C. Sagui, J. Board, T. Darden, Efficient particle-mesh ewald based approach to fixed and induced dipolar interactions, J. Chem. Phys. 113 (2000) 10913e10927. [42] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, Numerical-integration of cartesian equations of motion of a system with constraints e molecular-dynamics of nalkANES, J. Comput. Phys. 23 (1977) 327e341. [43] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys. 81 (1984) 3684e3690. [44] G.M. Torrie, J.P. Valleau, Nonphysical sampling distribution in Monte Carlo free-energy estimation: umbrella sampling, J. Comput. Phys. 23 (1977) 187e199. [45] J.G. Kirkwood, Statistical mechanics of fluid mixtures, J. Chem. Phys. 3 (1935) 300e313. [46] S. Kumar, D. Bouzida, R.H. Swendsen, P.A. Kollman, J.M. Rosenberg, The weighted histogram analysis method for free-energy calculations on biomolecules .1. The method, J. Comput. Chem. 13 (1992) 1011e1021. [47] B. Roux, The calculation of the potential of mean force using computer simulations, Comput. Phys. Commu. 91 (1995) 275e282. [48] B.C. Garrett, D.G. Truhlar, Generalized transition state theory. Quantum effects for collinear reactions of hydrogen molecules and isotopically substituted hydrogen molecules, J. Chem. Phys. 83 (1979) 1079e1112.
Contents lists available at ScienceDirect
Chemico-Biological Interactions journal homepage: www.elsevier.com/locate/chembioint
Free energy profiles of cocaine esterase-cocaine binding process by molecular dynamics and potential of mean force simulations Yuxin Zhang a, b, c, 1, Xiaoqin Huang b, 1, Keli Han a, Fang Zheng b, c, Chang-Guo Zhan b, c, * a
State Key Laboratory of Molecular Reaction Dynamics, Dalian Institute of Chemical Physics, Chinese Academy of Science, Zhongshan Road 457, Dalian 116023, PR China b Molecular Modeling and Biopharmaceutical Center, College of Pharmacy, University of Kentucky, 789 South Limestone Street, Lexington, KY 40536, USA c Department of Pharmaceutical Sciences, College of Pharmacy, University of Kentucky, 789 South Limestone Street, Lexington, KY 40536, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 9 January 2016 Received in revised form 29 April 2016 Accepted 5 May 2016 Available online 6 May 2016
The combined molecular dynamics (MD) and potential of mean force (PMF) simulations have been performed to determine the free energy profile of the CocE)-(þ)-cocaine binding process in comparison with that of the corresponding CocE-()-cocaine binding process. According to the MD simulations, the equilibrium CocE-(þ)-cocaine binding mode is similar to the CocE-()-cocaine binding mode. However, based on the simulated free energy profiles, a significant free energy barrier (~5 kcal/mol) exists in the CocE-(þ)-cocaine binding process whereas no obvious free energy barrier exists in the CocE-()-cocaine binding process, although the free energy barrier of ~5 kcal/mol is not high enough to really slow down the CocE-(þ)-cocaine binding process. In addition, the obtained free energy profiles also demonstrate that (þ)-cocaine and ()-cocaine have very close binding free energies with CocE, with a negligible difference (~0.2 kcal/mol), which is qualitatively consistent with the nearly same experimental KM values of the CocE enzyme for (þ)-cocaine and ()-cocaine. The consistency between the computational results and available experimental data suggests that the mechanistic insights obtained from this study are reasonable. © 2016 Elsevier Ireland Ltd. All rights reserved.
Keywords: Esterase Cocaine Binding process Energy barrier
1. Introduction Cocaine may have two enantiomers, ()-cocaine and (þ)-cocaine (see Fig. 1). The naturally occurring and biologically active ()-cocaine is known as the most addictive substance abused by millions of people worldwide [1e4]. Despite of decades of research efforts, there is no FDA-approved anti-cocaine therapeutic agent specific for cocaine addiction or overdose [4,5]. As well known, cocaine exerts its effects on central nervous system (CNS) through blocking the reuptake of neurotransmitter dopamine, thus potentiating the effects of dopamine in the synapse [6]. Traditional pharmacodynamic approaches have failed to produce a therapeutically useful small-molecule drug because of the difficulties inherent in blocking a blocker like cocaine. As an alternative,
* Corresponding author. Department of Pharmaceutical Sciences, College of Pharmacy, University of Kentucky, 789 South Limestone Street, Lexington, KY 40536, USA. E-mail address: [email protected] (C.-G. Zhan). 1 These authors contributed equally to this work. http://dx.doi.org/10.1016/j.cbi.2016.05.011 0009-2797/© 2016 Elsevier Ireland Ltd. All rights reserved.
pharmacokinetic approach using an efficient cocaine-metabolizing enzyme has been recognized as a promising strategy for therapeutic treatment of cocaine overdose and abuse. The enzyme strategy aims at accelerating cocaine metabolism (particularly through catalyzing cocaine benzoyl ester hydrolysis) and, therefore, eliminating cocaine rapidly in the plasma [1,7e9]. For this purpose, cocaine esterase (CocE) [10e13] is a promising choice as a potential anti-cocaine agent for therapeutic treatment of cocaine overdose and abuse, because CocE is the most efficient natural enzyme against ()-cocaine [12]. Native CocE is a bacterial enzyme with a catalytic center similar to that of butyrylcholinesterase (BChE) [14,15], but the catalytic efficiency of CocE [12] is two orders of magnitude higher than that of wild-type BChE against ()-cocaine [16]. Native CocE and the designed thermostable mutants are capable of protecting against cocaine-induced lethality [12,13,17,18]. Thus, both BChE and CocE are interesting protein scaffolds that may be used to design enzyme mutants with improved catalytic efficiency against ()-cocaine. To rationally redesign an enzyme (BChE or CocE) for development of anti-cocaine therapeutics, one first needs to understand the detailed catalytic mechanisms of the enzyme against both
Y. Zhang et al. / Chemico-Biological Interactions 259 (2016) 142e147
143
2. Methods 2.1. MD simulation
Fig. 1. Molecular structures of cocaine enantiomers.
()-cocaine and (þ)-cocaine and understand their similarities and differences in the enzyme-substrate binding and subsequent transformation processes. Previously reported studies [7,14,15,19,20] revealed that BChE- and CocE-catalyzed hydrolyses of ()-cocaine and (þ)-cocaine share a similar catalytic reaction pathway for the chemical reaction process. For each of the enzymatic reactions, the chemical reaction process always includes both the acylation and deacylation stages, regardless of the enzyme (CocE or BChE) and substrate (()-cocaine or (þ)-cocaine). There is also no significant difference in Michaelis-Menten constant (KM) between BChE-catalyzed ()-cocaine hydrolysis and BChEcatalyzed (þ)-cocaine hydrolysis [21]. However, the enzymesubstrate binding processes are remarkably different [22,23]. For example, it has been known that the enzyme-substrate binding process is the rate-determining step for BChE-catalyzed hydrolysis of ()-cocaine, whereas the chemical reaction process is the ratedetermining step for BChE-catalyzed hydrolysis of (þ)-cocaine [16,24,25]. Based on the understanding of the mechanistic difference, we were able to perform structure and mechanism based BChE redesign, and successfully designed and discovered new BChE mutants with considerably improved catalytic efficiency against ()-cocaine [25e32]. It has also been known that the CocE-()-cocaine binding process [33] is remarkably different from the BChE-()-cocaine binding process [23]. BChE-()-cocaine binding process is associated with a very high energy barrier [22,23], whereas the CocE()-cocaine binding process is barrier-less [33]. For the remaining mechanistic question to be addressed, it is unknown whether the CocE-(þ)-cocaine binding process is also different from the CocE()-cocaine binding process, because no computational study has been reported on the CocE-(þ)-cocaine binding process. Understanding the difference/similarity between the CocE-(þ)-cocaine and CocE-()-cocaine binding processes would be valuable for guiding further anti-cocaine enzyme redesign efforts. In the present study, we carried out molecular dynamics (MD) and potential of mean force (PMF) simulations to determine the free energy profile of the CocE-(þ)-cocaine binding process in comparison with the corresponding CocE-()-cocaine binding free energy profile. The combined MD and PMF simulations have revealed some interesting difference and similarity in the free energy profile between the CocE-(þ)-cocaine and CocE-()-cocaine binding processes.
The CocE structure used in this study is a thermally stable CocE mutant (T172R/G173Q) designed in our previous study [12]. It has been known that the T172R/G173Q mutations can significantly stabilize the CocE structure without changing the catalytic activity. The same T172R/G173Q mutant CocE was also used to study the CocE-()-cocaine binding process [33]. For convenience, below we simply refer CocE to the T172R/G173Q mutant CocE, unless indicated explicitly otherwise. Initial structure of the CocE-(þ)-cocaine binding complex was prepared based on the available X-ray crystal structure [34] of CocE (PDB code 3I2F for the T172R/G173Q mutant CocE) and the results of our previous molecular docking and MD simulations on CocE binding with ()-cocaine [14,33]. The initial CocE-(þ)-cocaine binding structure was energy-minimized by using the Sander module of Amber 12 program [35] via a combined use of the steepest descent/conjugate gradient algorithms, with a convergence criterion of 0.01 kcal mol1 Å1, and the non-bonded cutoff distance was set to 10.0 Å. MD simulation was performed by using the Sander module of the Amber 12 program package. The general procedure of the MD simulation was similar to that used in our previously reported other computational studies [12,13,33,36]. In particular, the molecular mechanics force field parameters and the partial charges of (þ)-cocaine atoms were adopted directly from those developed in our previous studies on other complexes [12,13,33,36]. Briefly, the partial charges of (þ)-cocaine atoms were calculated by using the restrained electrostatic potential-fitting (RESP) protocol implemented in the Antechamber module of the Amber 12 program, following the electrostatic potential (ESP) calculation at ab initio HF/6-31G* level using Gaussian 03 program [37]. For the energy minimization and MD simulation, the CocE(þ)-cocaine complex was solvated in an orthorhombic box of TIP3P water molecules [38] with a minimum solvent-wall distance of 10 Å. Sodium counter ions (Naþ) were added to neutralize the solvated system. The solvated system was gradually heated to 298.15 K by the weak-coupling method [39] and equilibrated for 400 ps. During the MD simulations, a 10.0 Å non-bonded interaction cutoff was used and the non-bonded list was updated every 25 steps. The motion for the mass center of the system was removed every 1000 steps. The particle-mesh Ewald (PME) method [40,41] was used to treat long-range electrostatic interactions. The lengths of covalent bonds involving hydrogen atoms were fixed with the SHAKE algorithm [42], enabling the use of a 2-fs time step to numerically integrate the equations of motion. Finally, the production MD was kept running for ~2.0 ns with a periodic boundary condition in the NTP ensemble at T ¼ 298.15 K with Berendsen temperature coupling, and at P ¼ 1 atm with isotropic moleculebased scaling [40,43]. 2.2. Potential of mean force (PMF) simulation In order to explore the free energy profiles for the process of (þ)-cocaine binding with CocE, the PMF simulation was carried out by using umbrella-sampling [44] MD simulation. The classic PMF definition [45] can be represented by a function of reaction coordinate as
uðcÞ ¼ RT ln〈rðcÞ〉 UðcÞ þ F
(1)
in which r(c) is the probability density along the reaction coordinate c, R is the gas constant, T is the simulation temperature, U(c) is the biasing potential applied in the umbrella-sampling MD simulation, and F is the normalization constant. According to this
144
Y. Zhang et al. / Chemico-Biological Interactions 259 (2016) 142e147
approach, the reaction coordinate is usually divided into different regions, i.e., windows, and each window is sampled separately. A biasing (umbrella) potential, i.e. U(c), is applied for each window in order to obtain nearly uniform sampling of the potential energy surface. In the present study, the reaction coordinate was defined as the distance from the mass center of the non-hydrogen atoms of (þ)-cocaine to the mass center of the non-hydrogen atoms on the side chains of residues H87, V121, and L146 of CocE. The total number of windows was 73, with a window size of 0.3 Å, covering the reaction coordinate values from ~13.385 Å to ~34.985 Å. The biasing force constant applied in different windows of umbrellasampling was 10.0 kcal/(molÅ2), as used in our previous PMF simulations on the CocE-()-cocaine binding process [33]. For each umbrella-sampling window, the initial complex structure was selected from the last snapshot of the PMF simulations of the previous window. The selected structure for each window was first equilibrated for 200 ps and then kept running for 800 ps for production sampling. The frequency for data collection was set to 1 fs, which was the same as that of the time step of the umbrellasampling MD. Previous computational study [33] on the CocE()-cocaine binding process demonstrated that this PMF protocol led to converged PMF free energy profiles, as further longer MD and PMF simulations did not significantly change the obtained free energy profiles. After the umbrella-sampling MD simulations for all windows were completed, the data collected from separate simulation windows were combined along the reaction coordinate. These data were then used to calculate the PMF for the whole binding process with the weighed histogram analysis method (WHAM) [46,47] using the code developed by Alan Grossfield (http://membrane.urmc.rochester.edu/Software/WHAM/WHAM. html). Notably, the equation for the biasing potential implemented in the Amber program for the PMF simulations is U(c) ¼ kAmber (cc0)2, whereas the corresponding equation implemented in the WHAM program for the PMF free energy analysis is U(c) ¼ (kWHAM/2) (cc0)2. So, the equations implemented in the two programs are inconsistent in the force constant. Clearly, kWHAM ¼ 2kAmber. In this study, we used kAmber ¼ 10 kcal/ mol and kWHAM ¼ 20 kcal/mol. It should be noted that such a difference between kWHAM and kAmber was not noted in the previous analysis [33] of the PMF data obtained for the CocE-()-cocaine binding. As a result, the previous analysis [33] of the PMF data obtained for the CocE-()-cocaine binding was performed inconsistently by using kAmber ¼ 10 kcal/mol and kWHAM ¼ 10 kcal/mol. In this way, the previous analysis using the WHAM program systematically under-estimated the PMF free energy profile (and the binding free energy) for the CocE-()-cocaine binding by a factor of 2, although the systematical under-estimation did not change any of the conclusions in the previous study [33]. All of the MD and umbrella-sampling MD simulations were performed on a supercomputer (a DELL Cluster with 388 nodes or 4816 processors) at the University of Kentucky’s Computer Center.
Fig. 2. Tracked internuclear distances in the MD-simulated CocE-(þ)-cocaine binding complex. (A) C32(COC)/S117 (black color) represents the distance between the carbonyl carbon on the benzoyl group of (þ)-cocaine and the hydroxyl oxygen of the S117 side chain; O33(COC)/Y44 (red color) refers to the distance between the carbonyl oxygen on the benzoyl group of (þ)-cocaine and the hydroxyl hydrogen of the Y44 side chain; O33(COC)/Y118 (blue color) represents the distance between the carbonyl oxygen on the benzoyl group of (þ)-cocaine and the backbone hydrogen of Y118. (B) S117/H287 (black color) and H287—D259 (red color) represent the distances related to hydrogen bonds within the catalytic triad S117-H287-D259. Y44/W166 (blue color) represents the hydrogen bond distance between the hydroxyl oxygen on the Y44 side chain and the side chain of W166. The block colored distances are not visible because they overlap with the blue colored distances. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3. Results and discussion 3.1. CocE-(þ)-cocaine binding structure According to the MD-simulated CocE-(þ)-cocaine binding structure, (þ)-cocaine binding site is located on the interface of the three domains of CocE (domains I to III). Depicted in Fig. 2 are important distances in the simulated CocE-(þ)-cocaine complex, and depicted in Fig. 3 is the structure in the last snapshot of the MD simulated CocE-(þ)-cocaine complex, revealing the protein environment of (þ)-cocaine in the complex. Figs. 2 and 3 clearly show the important hydrogen-bonding interactions that were very stable
Fig. 3. MD-simulated CocE-(þ)-cocaine binding structure. The enzyme is shown as golden ribbons, and (þ)-cocaine is in ball-and-stick style. Key residues of the enzyme within 5 Å around the (þ)-cocaine molecule are shown in stick style and colored by atom types. Hydrogen bonds are represented as dashed lines with labeled distances (in unit of Å). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Y. Zhang et al. / Chemico-Biological Interactions 259 (2016) 142e147
during the MD simulation process. In the simulated CocE(þ)-cocaine complex (Fig. 3), the benzoyl ester group of (þ)-cocaine is located in a sub-binding site composed of hydrophobic residues W151, W166, L169, L407, F408, F261, and P150 of CocE, packing in parallel with the aromatic side chain of W166, and in perpendicular with the aromatic side chain of F261. As shown in Fig. 3, S117, H287, and D259 residues form a stable catalytic triad with a stable hydrogen bonding network, and the carbonyl carbon of (þ)-cocaine benzoyl ester group is about ~3.5 Å away from the hydroxyl oxygen of S117 side chain. ~3.5 Å is an ideal distance for the hydroxyl oxygen to initiate the nucleophilic attack on the carbonyl carbon, according to our previous computational studies on other ester hydrolysis reactions catalyzed by a serine esterase (with a similar catalytic triad) like BChE [7,24]. In addition, the carbonyl oxygen on the benzoyl ester group of (þ)-cocaine interacts with the hydroxyl group of Y44 side chain through hydrogen-bonding interactions. As seen in Figs. 2 and 3, there is a strong hydrogen bond between the carbonyl oxygen of (þ)-cocaine benzoyl ester group and the hydroxyl group of Y44 side chain, with an O … H distance of ~1.9 Å. This hydrogen bond is expected to become stronger and stronger in going from the prereaction CocE-()-cocaine binding to the transition state and intermediate due to the increasing negative charge in the carbonyl oxygen atom during the acylation process. Notably, the CocE(þ)-cocaine binding mode simulated in this study is generally similar to the corresponding CocE-()-cocaine binding mode found in our previous study [33]. 3.2. Free energy profile of the CocE-(þ)-cocaine binding process In order to calculate the free energy of the CocE-(þ)-cocaine binding, the PMF simulations were carried out starting from the MD-simulated CocE-(þ)-cocaine complex structure. Based on the data collected from the umbrella-sampling MD simulations, the PMF for the CocE-(þ)-cocaine binding was determined. The distance between the mass center of the non-hydrogen atoms of (þ)-cocaine and the mass center of the non-hydrogen atoms on the side chains of residues H87, V121, and L146 of the enzyme was used as the reaction coordinate for the PMF calculations. The selection of the reaction coordinate was based on the structural features of our modeled CocE-(þ)-cocaine binding structures in the present study, because we found that the direction of reaction coordinate roughly went through the central point of the active site of CocE to reach the molecular surface of the enzyme. For a fair comparison between the free energy profile of the CocE-(þ)-cocaine binding process and that of the CocE-()-cocaine binding process, we also analyzed the previously obtained PMF data (based on the use of the same computational protocol described in the present study) for the CocE-()-cocaine binding process with kWHAM ¼ 20 kcal/mol (because the PMF simulations [33] were carried out using kAmber ¼ 10 kcal/mol). Depicted in Fig. 4 are the PMF free energy profiles obtained for the CocE-(þ)-cocaine and CocE-()-cocaine binding processes. Based on the PMF free energy profiles depicted in Fig. 4, there is a significant free energy barrier (~2 kcal/mol) when the reaction coordinate c ¼ ~17.6 Å (associated with a transition state) in the CocE-(þ)-cocaine binding process, whereas we cannot identify an obvious free energy barrier in the CocE-()-cocaine binding process. The free energy barrier for the CocE-(þ)-cocaine binding process is associated with (þ)-cocaine rotation in the binding pocket. Fig. 5(A) depicts the CocE-(þ)-cocaine complex structure when c ¼ ~17.8 Å (right before the transition state), and Fig. 5(B) shows the CocE-(þ)-cocaine complex structure when c ¼ ~17.3 Å (right after the transition state). Further, in terms of the binding kinetics, a free energy barrier of
145
Fig. 4. Simulated free energy profile of the CocE-(þ)-cocaine binding process (A) in comparison with that of the CocE-()-cocaine binding process (B). The reaction coordinate was defined as the distance between the mass center of the non-hydrogen atoms of (þ)/()-cocaine and the mass center of the non-hydrogen atoms on the side chains of residues H87, V121, and L146 of the enzyme.
~2 kcal/mol is associated with a rate constant as large as ~3 1011 s1 (faster than the diffusion) according to the classical rate constant theory [48], suggesting that the binding process is much faster than the chemical reaction process (with a much lower catalytic rate constant kcat [10,12]). So, the binding process associated with a free energy barrier of ~2 kcal/mol is still very fast and, hence, easy to reach the thermodynamic equilibration of the CocE(þ)-cocaine binding. So, the free energy barrier of ~2 kcal/mol is not high enough to really slow down the enzymatic reaction process (compared to the much slower chemical reaction process), and the CocE-substrate binding process is not rate-determining for both (þ)-cocaine and ()-cocaine. In this sense, we may also say that the CocE-(þ)-cocaine and CocE-()-cocaine binding processes are similar, in contrast to the remarkable difference between BChE(þ)-cocaine and BChE-()-cocaine binding processes. It has been known that the BChE-()-cocaine binding process is ratedetermining, whereas the BChE-(þ)-cocaine binding process is not [22e24]. In terms of the equilibrium CocE-substrate binding, the free energy profiles depicted in Fig. 4 indicate that the calculated CocE-substrate binding free energy is 12.2 kcal/mol for (þ)-cocaine and 12.4 kcal/mol for ()-cocaine. The difference (~0.2 kcal/mol) between (þ)-cocaine and ()-cocaine is negligible, which is qualitatively consistent with the nearly same experimental KM values of the CocE enzyme for the two substrates (KM ¼ ~15 mM for (þ)-cocaine and KM ¼ ~13 mM for ()-cocaine) [20]. 4. Conclusion The combined MD and PMF simulations have led to determination of the free energy profile of the CocE-(þ)-cocaine binding process in comparison with that of the CocE-()-cocaine binding
146
Y. Zhang et al. / Chemico-Biological Interactions 259 (2016) 142e147
DA041115) and R01 grants (R01 DA035552, R01 DA032910, R01 DA013930, and R01 DA025100). The authors also acknowledge the Computer Center at University of Kentucky for supercomputing time on a Dell Supercomputer Cluster consisting of 388 nodes or 4816 processors. Transparency document Transparency document related to this article can be found online at http://dx.doi.org/10.1016/j.cbi.2016.05.011. References
Fig. 5. Simulated CocE-(þ)-cocaine binding structures right before and after the rotational transition state (when c ¼ ~17.6 Å): (A) the CocE-(þ)-cocaine binding structure when c ¼ ~17.8 Å (right before the transition state); (B) the CocE-(þ)-cocaine binding structure c ¼ ~17.3 Å (right after the transition state).
process. According to the obtained free energy profiles, a significant free energy barrier (~5 kcal/mol) exists in the CocE-(þ)-cocaine binding process, whereas no obvious free energy barrier exists in the CocE-()-cocaine binding process. The free energy barrier for the CocE-(þ)-cocaine binding process is associated with (þ)-cocaine rotation in the binding pocket. On the other hand, the free energy barrier of ~2 kcal/mol is not high enough to really slow down the CocE-(þ)-cocaine binding process, and the CocEsubstrate binding process should not be rate-determining for both (þ)-cocaine and ()-cocaine. The obtained free energy profiles also reveal that (þ)-cocaine and ()-cocaine have very close binding free energies with CocE. The negligible binding free energy difference (~0.2 kcal/mol) between (þ)-cocaine and ()-cocaine is qualitatively consistent with the nearly same experimental KM values of the CocE enzyme for (þ)-cocaine and ()-cocaine. The consistency between the computational and experimental data suggests that the mechanistic insights obtained from the computational simulations in this study are reasonable. Competing financial interests statement The authors declare that there is no conflict of interest for this work. Acknowledgments This work was supported by the National Institutes of Health (NIH) through the NIDA Translational Avant-Garde Award (UH2
[1] D.W. Landry, K. Zhao, G.X. Yang, M. Glickman, T.M. Georgiadis, Antibodycatalyzed degradation of cocaine, Science 259 (1993) 1899e1901. [2] A.L. Milton, B.J. Everitt, Wiping drug memories, Science 336 (2012) 167e168. [3] K.D. Ersche, P.S. Jones, G.B. Williams, A.J. Turton, T.W. Robbins, E.T. Bullmore, Abnormal brain structure implicated in stimulant drug addiction, Science 335 (2012) 601e604. [4] F. Zheng, C.-G. Zhan, Are pharmacokinetic approaches feasible for treatment of cocaine addiction and overdose? Future Med. Chem. 4 (2012) 125e128. [5] D.A. Gorelick, Pharmacokinetic approaches to treatment of drug addiction, Exp. Rev. Clin. Pharmacol. 1 (2008) 277e290. [6] R. Chen, M.R. Tilley, H. Wei, F. Zhou, F.-M. Zhou, S. Ching, N. Quan, R.L. Stephens, E.R. Hill, T. Nottoli, D.D. Han, H.H. Gu, Abolished cocaine reward in mice with a cocaine-insensitive dopamine transporter, Proc. Natl. Acad. Sci. U. S. A. 103 (2006) 9333e9338. [7] C.-G. Zhan, F. Zheng, D.W. Landry, Fundamental reaction mechanism for cocaine hydrolysis in human butyrylcholinesterase, J. Am. Chem. Soc. 125 (2003) 2462e2474. [8] F. Zheng, C.-G. Zhan, Modeling of pharmacokinetics of cocaine in human reveals the feasibility for development of enzyme therapies for drugs of abuse, PLoS Comput. Biol. 8 (2012) e1002610. [9] X. Chen, L. Xue, S. Hou, Z. Jin, T. Zhang, F. Zheng, C.-G. Zhan, Long-acting cocaine hydrolase for addiction therapy, Proc. Natl. Acad. Sci. U. S. A. 113 (2016) 422e427. [10] N.A. Larsen, J.M. Turner, J. Stevens, S.J. Rosser, A. Basran, R.A. Lerner, N.C. Bruce, I.A. Wilson, Crystal structure of a bacterial cocaine esterase, Nat. Struct. Mol. Biol. 9 (2002) 17e21. [11] M.C. Ko, L.D. Bowen, D. Narasimhan, A.A. Berlin, N.W. Lukacs, R.K. Sunahara, Z.D. Cooper, J.H. Woods, Cocaine esterase: interactions with cocaine and immune responses in mice, J. Pharmacol. Exp. Ther. 320 (2007) 926e933. [12] D. Gao, D.L. Narasimhan, J. Macdonald, M.-C. Ko, D.W. Landry, J.H. Woods, R.K. Sunahara, C.-G. Zhan, Thermostable variants of cocaine esterase for longtime protection against cocaine toxicity, Mol. Pharmacol. 75 (2009) 318e323. [13] L. Fang, K.M. Chow, S. Hou, L. Xue, D.W. Rodgers, F. Zheng, C.-G. Zhan, Rational design, preparation, and characterization of a therapeutic enzyme mutant with improved stability and function for cocaine detoxification, ACS Chem. Biol. 9 (2014) 1764e1772. [14] J.J. Liu, A. Hamza, C.-G. Zhan, Fundamental reaction mechanism free energy profile for (-)-cocaine hydrolysis catalyzed by cocaine esterase, J. Am. Chem. Soc. 131 (2009) 11964e11975. [15] J. Liu, C.-G. Zhan, Reaction pathway and free energy profile for cocaine hydrolase-catalyzed hydrolysis of (e)-cocaine, J. Chem. Theory Comput. 8 (2012) 1426e1435. [16] H. Sun, J. El Yazal, O. Lockridge, L.M. Schopfer, S. Brimijoin, Y.-P. Pang, Predicted michaelis-menten complexes of cocaine-butyrylcholinesterase: engineering effective butyrylcholinesterase mutants for cocaine detoxication, J. Biol. Chem. 276 (2001) 9330e9336. [17] G.T. Collins, R.L. Brim, D. Narasimhan, M.-C. Ko, R.K. Sunahara, C.-G. Zhan, J.H. Woods, Cocaine esterase prevents cocaine-induced toxicity and the ongoing intravenous self-administration of cocaine in rats, J. Pharmacol. Exp. Ther. 331 (2009) 445e455. [18] R.L. Brim, M.R. Nance, D.W. Youngstrom, D. Narasimhan, C.-G. Zhan, J.J. Tesmer, R.K. Sunahara, J.H. Woods, A thermally stable form of bacterial cocaine esterase: a potential therapeutic agent for treatment of cocaine abuse, Mol. Pharmacol. 77 (2010) 593e600. [19] C.-G. Zhan, D. Gao, Catalytic mechanism and energy barriers for butyrylcholinesterase-catalyzed hydrolysis of cocaine, Biophys. J. 89 (2005) 3863e3872. [20] J.J. Liu, X.Y. Zhao, W.C. Yang, C.-G. Zhan, Reaction mechanism for cocaine esterase-catalyzed hydrolyses of (þ)- and (-)-cocaine: unexpected common rate-determining step, J. Phys. Chem. B 115 (2011) 5017e5025. [21] L. Xue, S. Hou, W. Yang, L. Fang, F. Zheng, C.-G. Zhan, Catalytic activities of a cocaine hydrolase engineered from human butyrylcholinesterase against (þ)and (-)-cocaine, Chem. Biol. Interact. 203 (2013) 57e62. [22] X. Huang, Y. Pan, F. Zheng, C.-G. Zhan, Reaction pathway and free energy profile for pre-chemical reaction step of human butyrylcholinesterasecatalyzed hydrolysis of (-)-cocaine by combined targeted molecular dynamics and potential of mean force simulations, J. Phys. Chem. B 114 (2010) 13545e13554.
Y. Zhang et al. / Chemico-Biological Interactions 259 (2016) 142e147 [23] X. Huang, F. Zheng, C.-G. Zhan, Human butyrylcholinesterase-cocaine binding pathway and free energy profiles by molecular dynamics and potential of mean force simulations, J. Phys. Chem. B 115 (2011) 11254e11260. [24] A. Hamza, H. Cho, H.-H. Tai, C.-G. Zhan, Molecular dynamics simulation of cocaine binding with human butyrylcholinesterase and its mutants, J. Phys. Chem. B 109 (2005) 4776e4782. [25] Y. Pan, D. Gao, W. Yang, H. Cho, G. Yang, H.-H. Tai, C.G. Zhan, Computational redesign of human butyrylcholinesterase for anticocaine medication, Proc. Natl. Acad. Sci. U. S. A. 102 (2005) 16656e16661. [26] F. Zheng, W.C. Yang, M.C. Ko, J.J. Liu, H. Cho, D.Q. Gao, M. Tong, H.H. Tai, J.H. Woods, C.-G. Zhan, Most efficient cocaine hydrolase designed by virtual screening of transition states, J. Am. Chem. Soc. 130 (2008) 12148e12155. [27] F. Zheng, W. Yang, L. Xue, S. Hou, J. Liu, C.-G. Zhan, Design of high-activity mutants of human butyrylcholinesterase against (-)-cocaine: structural and energetic factors affecting the catalytic efficiency, Biochemistry 49 (2010) 9113e9119. [28] L. Xue, M.-C. Ko, M. Tong, W. Yang, S. Hou, L. Fang, J. Liu, F. Zheng, J.H. Woods, H.-H. Tai, C.-G. Zhan, Design, preparation, and characterization of high-activity mutants of human butyrylcholinesterase specific for detoxification of cocaine, Mol. Pharmacol. 79 (2011) 290e297. [29] W. Yang, Y. Pan, L. Fang, D. Gao, F. Zheng, C.-G. Zhan, Free-energy perturbation simulation on transition states and high-activity mutants of human butyrylcholinesterase for (-)-cocaine hydrolysis, J. Phys. Chem. B 114 (2010) 10889e10896. [30] W. Yang, Y. Pan, F. Zheng, H. Cho, H.-H. Tai, C.-G. Zhan, Free-energy perturbation simulation on transition states and redesign of butyrylcholinesterase, Biophys. J. 96 (2009) 1931e1938. [31] F. Zheng, L. Xue, S. Hou, J. Liu, M. Zhan, W. Yang, C.-G. Zhan, A highly efficient cocaine-detoxifying enzyme obtained by computational design, Nat. Commun. 5 (2014) 3457, http://dx.doi.org/10.1388/ncomms4457. [32] D. Gao, H. Cho, W. Yang, Y. Pan, G. Yang, H.H. Tai, C.-G. Zhan, Computational design of a human butyrylcholinesterase mutant for accelerating cocaine hydrolysis based on the transition-state simulation, Angew. Chem. Int. Ed. 45 (2006) 653e657. [33] X. Huang, X. Zhao, F. Zheng, C.-G. Zhan, Cocaine esterase-cocaine binding process and the free energy profiles by molecular dynamics and potential of mean force simulations, J. Phys. Chem. B 116 (2012) 3361e3368. [34] D. Narasimhan, M.R. Nance, D. Gao, M.-C. Ko, J. Macdonald, P. Tamburi, D. Yoon, D.M. Landry, J.H. Woods, C.-G. Zhan, J.J.G. Tesmer, R.K. Sunahara, Structural analysis of thermostabilizing mutations of cocaine esterase, Protein. Eng. Des. Sel. 23 (2010) 537e547. [35] D.A. Case, T.A. Darden, T.E. Cheatham Iii, C.L. Simmerling, J. Wang, R.E. Duke, R. Luo, R.C. Walker, W. Zhang, K.M. Merz, B. Roberts, S. Hayik, A. Roitberg, ry, K.F. Wong, F. Paesani, J. Vanicek, G. Seabra, J. Swails, A.W. Goetz, I. Kolossva R.M. Wolf, J. Liu, X. Wu, S.R. Brozell, T. Steinbrecher, H. Gohlke, Q. Cai, X. Ye, J. Wang, M.J. Hsieh, G. Cui, D.R. Roe, D.H. Mathews, M.G. Seetin, R. SalomonFerrer, C. Sagui, V. Babin, T. Luchko, S. Gusarov, A. Kovalenko, P.A. Kollman, AMBER 12, University of California, San Francisco, 2012.
147
[36] J. Liu, A. Hamza, C.-G. Zhan, Fundamental reaction mechanism free energy profile for ()-cocaine hydrolysis catalyzed by cocaine esterase, J. Am. Chem. Soc. 131 (2009) 11964e11975. [37] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, J.A. Montgomery Jr., T. Vreven, K.N. Kudin, J.C. Burant, J.M. Millam, S.S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G.A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J.E. Knox, H.P. Hratchian, J.B. Cross, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, P.Y. Ayala, K. Morokuma, G.A. Voth, P. Salvador, J.J. Dannenberg, V.G. Zakrzewski, S. Dapprich, A.D. Daniels, M.C. Strain, O. Farkas, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J.V. Ortiz, Q. Cui, A.G. Baboul, S. Clifford, J. Cioslowski, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, M.W. Wong, C. Gonzalez, J.A. Pople, Gaussian 03, Rev A.1, Gaussian Inc., Pittsburgh, PA, 2003. [38] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys. 79 (1983) 926e935. [39] T. Morishita, Fluctuation formulas in molecular dynamics simulations with the weak coupling heat bath, J. Chem. Phys. 113 (2000) 2976e2982. [40] T. Darden, D. York, L. Pedersen, Particle mesh ewalddan N_log(N) method for ewald sums in large systems, J. Chem. Phys. 98 (1993) 10089e10092. [41] A. Toukmaji, C. Sagui, J. Board, T. Darden, Efficient particle-mesh ewald based approach to fixed and induced dipolar interactions, J. Chem. Phys. 113 (2000) 10913e10927. [42] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, Numerical-integration of cartesian equations of motion of a system with constraints e molecular-dynamics of nalkANES, J. Comput. Phys. 23 (1977) 327e341. [43] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys. 81 (1984) 3684e3690. [44] G.M. Torrie, J.P. Valleau, Nonphysical sampling distribution in Monte Carlo free-energy estimation: umbrella sampling, J. Comput. Phys. 23 (1977) 187e199. [45] J.G. Kirkwood, Statistical mechanics of fluid mixtures, J. Chem. Phys. 3 (1935) 300e313. [46] S. Kumar, D. Bouzida, R.H. Swendsen, P.A. Kollman, J.M. Rosenberg, The weighted histogram analysis method for free-energy calculations on biomolecules .1. The method, J. Comput. Chem. 13 (1992) 1011e1021. [47] B. Roux, The calculation of the potential of mean force using computer simulations, Comput. Phys. Commu. 91 (1995) 275e282. [48] B.C. Garrett, D.G. Truhlar, Generalized transition state theory. Quantum effects for collinear reactions of hydrogen molecules and isotopically substituted hydrogen molecules, J. Chem. Phys. 83 (1979) 1079e1112.