Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Bioenergetics. Author manuscript; available in PMC 2017 September 21.
Published in final edited form as:
Published online 2017 January 8. doi:  10.4172/2167-7662.1000145
PMCID: PMC5608097

Ensemble Molecular Dynamics of a Protein-Ligand Complex: Residual Inhibitor Entropy Enhances Drug Potency in Butyrylcholinesterase


Butyrylcholinesterase is a key enzyme that catalyzes the hydrolysis of the neurotransmitter acetylcholine and shows an increased activity in patients suffering from Alzheimer’s disease (AD), making this enzyme a primary target in treating AD. Central to this problem, and to similar scenarios involving biomolecular recognition, is our understanding of the nature of the protein-ligand complex. The butyrylcholinesterase enzyme was studied via all-atom, explicit solvent, ensemble molecular dynamics simulations sans inhibitor and in the presence of three dialkyl phenyl phosphate inhibitors of known potency to a cumulative sampling of over 40 μs. Following the relaxation of these ensembles to conformational equilibria, binding modes for each inhibitor were identified. While classical models, which assume significant reduction in protein and ligand conformational entropies, continue to be favored in contemporary studies, our observations contradict those assumptions: bound ligands occupy many conformational states, thereby stabilizing the complex, while also promoting protein flexibility.

Keywords: BChE, Distributed computing, Phenyl phosphate, Molecular simulation, Docking


Acetylcholinesterase (AChE) is highly specific to the neurotransmitter acetylcholine, which the enzyme hydrolyzes, resulting in the termination of synaptic neurotransmission [1,2]. In contrast, butyrylcholinesterase (BChE), shown in (Figure 1), is a promiscuous enzyme that can act as an AChE-substitute in vivo and can hydrolyze various cholines, acyl cholines, acyl thiocholines, succinyl cholines, organophosphates, and acetanilides [36]. These enzymes, which play physiochemically distinct roles both in neurotransmission and in cellular differentiation and development [7], have thus been targeted as biosensors and bioscavengers that can detect and detoxify a myriad of organic poisons and pesticides [8,3,9].

Figure 1
X-ray structure of BChE (1P0I.pdb) with the binding pocket magnified and a schematic of the DAPP inhibitors studied. Active site residues are color coded for easy identification including the catalytic triad (yellow), the oxyanion hole (orange), the choline ...

They have also been targeted in treating a number of human health conditions including glaucoma, myasthenia gravis, and various central nervous system disorders such as traumatic brain injury, Down syndrome, and Parkinson dementia [1012]. Often noted are the roles that these enzymes play in Alzheimer’s disease (AD), with AChE concentrations decreasing as the disease progresses, and BChE levels increasing to take up the role of hydrolysis in cholinergic neurons, concomitant with increasing quantities of amyloid-rich neural plaques and tangles [13].

The development of natural and synthetic cholinesterase inhibitors is thus a rapidly evolving field, and understanding both the physical interactions upon which protein-ligand binding is dependent, and the dynamics inherent to such events, is central to future progress in all areas of biomolecular recognition.

A primary supposition of classical models of protein-ligand binding is that the binding event results in a single, lowest-energy protein-ligand configuration [14], which implicitly assumes that both protein and bound ligand suffer significant penalties in configurational entropy upon binding, thereby requiring an enthalpic counterbalance to overcome this entropy loss. In contrast, we observe a significant number of inhibitor-accessible binding modes that not only endow the inhibitor with unpredicted, residual conformational entropy, but also allow for continued protein flexibility after ligand binding.

We report herein our computational study of the inhibition of BchE, which is targeted for moderate to severe symptoms of AD [2,13]. As in AchE, the Ser-His-Glu catalytic triad of BchE (yellow in (Figure 1) inset) is located near the bottom of a binding pocket pprox. 20 Å deep, across which exists an electrostatic gradient. This active site gorge is lined with a number of aromatic and aliphatic residues, and previous studies have characterized the importance of a number of chemical groups within this pocket including the oxyanion hole, the acyl and choline binding sites, and the peripheral anionic site [1518] (orange, blue, green, and red, respectively, in the (Figure 1) inset). We employ three chemically similar, reversible dialkyl phenyl phosphate (DAPP) inhibitors, PO4(Ph)R2 (where R=methyl for DAPP1, n-propyl for DAPP3, and n-pentyl for DAPP5). These inhibitors have been characterized experimentally, displaying binding affinities that increase with alkyl chain length [15], and which form the base units of intriguing dimeric inhibitors that show improved potency [19].


All-atom molecular dynamics simulations of native BchE [16] sans inhibitor, and of the protein in complex with each of these DAPP inhibitors, were performed using GROMACS 3.3 [20] inside the Folding@Home distributed computing architecture [21]. The ~60 kDa protein was modeled using the AMBER-03 force field of Duan, et al. [22]. Inhibitors were modeled using the general AMBER force field (GAFF) with RESP charges derived at the 6-31G* level using RED Server [23]. After initial docking of each inhibitor, sodium ions were randomly placed in a cubic 100 Å periodic box centered on the protein to establish electro neutrality. Solvation of this system with ~30,200 TIP3P [24] explicit water molecules resulted in a total system size of nearly 100,000 atoms. Following annealing of the ionic solvent with the protein held fixed, 1000 simulations of each system were initiated. All simulations were performed in the NPT ensemble [25] at 1.0 atm and 300 K with switched cutoffs applied to van der Waals interactions between 8.0 and 10.0 Å, and electrostatic interactions beyond 12.0 Å treated via reaction field. A 2.0 fs timestep was used, with bonds involving hydrogen atoms constrained using the LINCS algorithm [26], and conformations stored every 100 ps. We stress that no artificial or biasing potentials or restraints were applied to any portion of these simulated systems.

A total sampling time of ~7 μs per BchE-DAPP complex and ~20 μs for the enzyme sans inhibitor was collected, yielding a cumulative sampling of over 40 μs. In each simulated ensemble, the size and native structure of the protein was completely maintained, as demonstrated in (Table 1). For each DAPP inhibitor, the P-O-Ph oxygen was taken as the center-of-geometry (COG). All configurations in the resulting data sets were then analyzed by first aligning the protein to the initial (reference) structure and then characterizing the inhibitor position and orientation using a 15-dimensional vector composed of (a) the vector defining the inhibitor COG position relative to the protein center-of-mass, (b) a vector defining the directional axis through the COG and phenyl para-carbon, (c) the vector normal to the phenyl ring plane, and (d) vectors defining the directions of the alkyl groups relative to the COG. For each conformation, this 15D vector was used as the basis for K-means clustering [27] to identify inhibitor binding modes.

Table 1
Ensemble averaged structural metrics for BChE alone and in BChE-DAPP complexes.

Results and Discussion

Based on the population of each binding mode monitored over time, as illustrated for DAPP5 in (Figure 2), conformational equilibrium was determined to occur at or prior to 6.0 ns in each ensemble. Following K-means clustering, the transition matrices for all DAPP ensemble, representing moves from each binding mode i to each binding mode j after this equilibration period, were found to demonstrate both time-independence and detailed balance. Thermodynamic quantities were thus evaluated using only data beyond this point. (Table 2) shows the number of binding modes (Nbind), the percentage of equilibrium in which docked inhibitor conformations were observed (%bind), and the configurational entropy associated with each bound inhibitor (TSbind), which was calculated using the statistical weight of each observed binding mode, Sbind =−R Σ w(i) ln w(i), as defined by Chandler [29].

Figure 2
Populations of observed DAPP5 binding modes versus time, with binding mode populations split into three panels for visual clarity.
Table 2
DAPP inhibitor binding properties.

Also shown in (Table 2) are the experimentally observed inhibitor dissociation constants (KI) for these three inhibitors [15], which are the equilibrium constants for the undocking process BChE·I → BChE+I, and which we can thus approximate based on the number of unbinding events observed in our ensemble simulations as KI=[BChE] [I]/[BChE·I].

Significant dissociation of the DAPP1 inhibitor was observed, with the inhibitor re-entering the pocket in some simulations and leaving the enzyme entirely or interacting with the protein surface in others, making simple approximation of the DAPP1 dissociation constant intractable.

For DAPP3 and DAPP5, however, we observed 9 and 1 unbinding events, respectively, translating to approximate KI values of 0.082 and 0.001. These approximate values are well in-line with the experimental values tabulated above and validate both our simulation methodology and the inherent binding strength of these DAPPs in silico.

As shown in (Table 2), the three DAPP inhibitors studied herein sample numerous binding modes, which enhances the conformational entropy of the inhibitor after binding. This residual inhibitor entropy, TSbind, which is not accounted for in classical models, increases from ~2 kT to ~3 kT going from DAPP1 to DAPP5, increasing linearly with the log of the binding strength.

Moreover, while classical models assume that the configurational entropy of the protein will also decrease upon binding, we observe no sign of decreasing protein flexibility.

Indeed, as illustrated in (Figure 3), root-mean-squared fluctuations (RMSF) of the protein, per residue, show no significant change in protein flexibility within the binding pocket or elsewhere, and suggests quite the opposite: the RMSF of some residues in and near the binding pocket increase slightly in the presence of the more potent DAPP3 and DAPP5 inhibitors.

Figure 3
Top: The root-mean-squared fluctuation (RMSF) for each residue is shown for the enzyme sans inhibitor (black), and for the enzyme in the presence of the DAPP1 (blue), DAPP3 (green), and DAPP5 (red) inhibitors. Bottom: The percentage change in RMSF per ...

We conclude that the diverse binding modes observed herein not only enhance inhibitor entropy, but also promote protein flexibility and, by extension, protein configurational entropy. Structural representations of the 24 binding modes observed for the strong DAPP5 inhibitor, along with their relative binding free energies, are shown in (Figure 4).

Figure 4
Average binding conformations for the DAPP5 inhibitor, which is shown in stick mode with phenyl and alkyl carbons shown in cyan and phosphorus, oxygen, and hydrogen atoms shown in yellow, red, and white, respectively. Water molecules have been removed ...

These observations force us to question the applicability of classically-based docking models, which is well-supported by recent studies, experimental and computational alike, that have identified unexpected entropic contributions to a number of phenomena involving molecular recognition. For example, Mao et al. [30] applied scanning tunneling microscopy to study the binding of thioflavin T peptide to a prion peptide, identifying four binding modes of varying statistical weight and suggesting that more modes were possible.

Cramer and co-workers combined all-atom simulation and nuclear magnetic resonance measurements to study ubiquitin in an aqueous solution of free ligands, reporting that bound ligands accessed a number of favorable conformations, and suggesting only a moderate loss of entropy upon binding [31]. And Lee and co-workers applied thermochemical measurements and analysis of crystallographic data to examine the inhibition of HIV-1 protease, noting a degeneracy of inhibitor binding states that is enhanced via solvent anchoring of the inhibitor to the active site [32].

These studies, alongside a number of additional observations put forth in the last decade, have strongly emphasized the importance of conformational flexibility and entropy in protein-ligand complex formation and stabilization [3335], which are both accounted for in our all-atom ensemble simulations.

Moreover, a notable review by Mobley and Dill on the physics of ligand binding emphasized the notion that small changes in conformation can lead to large changes in binding affinity [36], which is well illustrated by our ensemble simulations. (Figure 4) provides concise descriptions of the interactions between DAPP5 and the BChE active site gorge for each observed binding mode, which includes residues in nearly every “hot spot” within the active site gorge identified by Butini et al. using bioinformatical techniques [37].

While statistical treatments provide an attractive model by which to discretize ligand binding modes, thereby allowing the tabulation of specific conformational preferences, biomolecular recognition (a.k.a. the “docking problem”) is more accurately modeled as the diffusive sampling of a continuous and rugged free energy landscape; the protein-ligand complex is a fluid body that is driven between local free energy minima by thermal fluctuations and solvent interactions, and the average binding mode structures depicted in Figure 4 thus represent only the most populated regions of this continuous free energy landscape.


The residual inhibitor entropy provided by diverse binding conformations, as described here for DAPP inhibition of BChE, increases the stability of the protein-inhibitor complex beyond classically-derived models, while also promoting protein flexibility. We postulate that this observation is generalizable to all flexible ligand-receptor pairs, particularly those in which a large, chemically-rich binding site is available to the ligand. While classical models are useful in providing an elementary understanding of protein-ligand binding, an ever-growing number of observations have demonstrated that even a qualitatively-accurate description of the protein-ligand complex must include conformational flexibility and entropic contributions, both of which must be accounted for quantitatively if future studies involving biomolecular recognition, docking, and drug design are to be successful.

Still, a number of questions remain, and a natural step forward in our analysis is a rigorous statistical treatment of the interactions described in (Figure 4), as well as an assessment of the role of water in these binding interactions, which has become a prominent factor in discussions of the physics of protein-ligand binding in recent years [38,39]. Additional future directions include an evaluation of the role of inhibitor chemistry in defining binding mode diversity via ensemble simulations of a large and chemically-disparate set of inhibitors, and characterization of the mechanism of protein-ligand association, which has once again become a pronounced point of discussion in the biochemical and biophysical communities [14,36].


The authors thank the Folding@Home volunteers worldwide who contributed invaluable processor time to this effort. This work was supported by a Research Corporation Cottrell College Science Award (#7840) and the MARC U*STAR program (AR), funded by NIH/NIGMS #T34GM008074. SC thanks Women & Philanthropy for providing undergraduate research scholarships to support her efforts.


Associated Content

Selected simulation movies are available at


1. Lleó A, Greenberg SM, Growdon JH. Current pharmacotherapy for Alzheimer’s disease. Annu Rev Med. 2006;57:513–533. [PubMed]
2. Darvesh S, Hopkins DA, Geula C. Neurobiology of butyrylcholinesterase. Nat Rev Neurosci. 2003;4:131–138. [PubMed]
3. Masson P, Lockridge O. Butyrylcholinesterase for protection from organophosphorus poisons: catalytic complexities and hysteretic behavior. Arch Biochem Biophys. 2010;494:107–120. [PMC free article] [PubMed]
4. Chen X, Fang L, Liu JJ, Zhan CG. Reaction Pathway and Free Energy Profiles for Butyrylcholinesterase-Catalyzed Hydrolysis of Acetylthiocholine. Biochemistry. 2012;51:1297–1305. [PMC free article] [PubMed]
5. Huang XQ, Zheng F, Zhan CG. Human Butyrylcholinesterase-Cocaine Binding Pathway and Free Energy Profiles by Molecular Dynamics and Potential of Mean Force Simulations. J Phys Chem B. 2011;115:11254–11260. [PMC free article] [PubMed]
6. Fang L, Pan YM, Muzyka JL, Zhan CG. Active Site Gating and Substrate Specificity of Butyrylcholinesterase and Acetylcholinesterase: Insights from Molecular Dynamics Simulations. J Phys Chem B. 2011;115(27):8797–8805. [PMC free article] [PubMed]
7. Sperling LE, Steinert G, Boutter J, Landgraf D, Hescheler J, et al. Characterisation of cholinesterase expression during murine embryonic stem cell differentiation. Chem-Biol Interact. 2008;175:156–160. [PubMed]
8. Pohanka M, Musilek K, Kuca K. Progress of biosensors based on cholinesterase inhibition. Curr Med Chem. 2009;16:1790–1798. [PubMed]
9. Van Dyk JS, Pletschke B. Review on the use of enzymes for the detection of organochlorine, organophosphate and carbamate pesticides in the environment. Chemosphere. 2011;82:291–307. [PubMed]
10. Lindner A, Schalke B, Toyka KV. Outcome in juvenile-onset myasthenia gravis: a retrospective study with long-term follow-up of 79 patients. J Neurol. 1997;244:515–520. [PubMed]
11. Orhan IE, Orhan G, Gurkas E. An overview on natural cholinesterase inhibitors--a multi-targeted drug class--and their mass production. Mini Rev Med Chem. 2011;11:836–842. [PubMed]
12. Giacobini E. Cholinesterase inhibitors: new roles and therapeutic alternatives. Pharmacol Res. 2004;50:433–440. [PubMed]
13. Greig NH, Utsuki T, Ingram DK, Wang Y, Pepeu G, et al. Selective butyrylcholinesterase inhibition elevates brain acetylcholine, augments learning and lowers Alzheimer beta-amyloid peptide in rodent. Proc Natl Acad Sci U S A. 2005;102:17213–17218. [PubMed]
14. Vértessy BG, Orosz F. From “fluctuation fit” to “conformational selection”: evolution, rediscovery, and integration of a concept. Bioessays. 2011;33:30–34. [PubMed]
15. Law K-S, Acey RA, Smith CR, Benton DA, Soroushian S, et al. Dialkyl phenyl phosphates as novel selective inhibitors of butyrylcholinesterase. Biochem Biophys Res Commun. 2007;355:371–378. [PubMed]
16. Nicolet Y, Lockridge O, Patrick Masson P, Fontecilla-Camps JC, Nachon F. Crystal Structure of Human Butyrylcholinesterase and of Its Complexes with Substrate and Products. J Biol Chem. 2003;278:41141–41147. [PubMed]
17. Dong Y, Ming Zha X. Synthesis, Biological Evaluation and Molecular Modeling of (E)- 3-Propargylene-, 3-Dihydro-2H-Indol-2-Ones as Acetyl- and Butyrylcholinesterase Inhibitors for the Treatment of Alzheimer’s Disease. Medicinal chemistry (Los Angeles) 2016;06(06)
18. Jin H. Acetylcholinesterase and Butyrylcholinesterase Inhibitory Properties of Functionalized Tetrahydroacridines and Related Analogs. Medicinal chemistry (Los Angeles) 2014;4(10)
19. Stein JM. The effect of adrenaline and of alpha- and beta-adrenergic blocking agents on ATP concentration and on incorporation of 32Pi into ATP in rat fat cells. Biochem Pharmacol. 1975;24:1659–1662. [PubMed]
20. Lindahl E, Hess B, van der Spoel D. GROMACS 3.0: a package for molecular simulation and trajectory analysis. J Mol Model. 2001;7:306–317.
21. Zagrovic B, Sorin EJ, Pande V. Beta-hairpin folding simulations in atomistic detail using an implicit solvent model. J Mol Biol. 2001;313:151–169. [PubMed]
22. Duan Y, Wu C, Chowdhury S, Lee MC, Xiong G, et al. A Point-Charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. J Comput Chem. 2003;24:1999–2012. [PubMed]
23. Vanquelef E, Simon S, Marquant G, Garcia E, Klimerak G, et al. RED Server: a web service for deriving RESP and ESP charges and building force field libraries for new molecules and molecular fragments. Nucleic Acids Res. 2011;39:W511–W517. [PMC free article] [PubMed]
24. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of simple potential functions for simulating liquid water. J Chem Phys. 1983;79:926–935.
25. Berendsen HJC, Postma JPM, Van Gunsteren WF, Dinola A, Haak J. Molecular-dynamics with coupling to an external bath. J Chem Phys. 1984;81:3684–3690.
26. Hess B, Bekker H, Berendsen HJC, Fraaije JGEM. LINCS: A Linear Constraint Solver for molecular simulations. J Comput Chem. 1997;18:1463–1472.
27. Sorin EJ, Pande VS. Exploring the helix-coil transition via all-atom equilibrium ensemble simulations. Biophys J. 2005;88:2472–2493. [PubMed]
28. Kabsch W, Sander C. Dictionary of protein secondary structure: pattern recognition of hydrogen-bonded and geometrical features. Biopolymers. 1983;22:2577–2637. [PubMed]
29. Chandler D. Statistical Mechanics Introduction to Modern Statistical Mechanics. Oxford University Press; New York: 1987. pp. 54–85.
30. Mao XB, Guo YY, Wang CX, Zhang M, Ma XJ, et al. Binding Modes of Thioflavin T Molecules to Prion Peptide Assemblies Identified by Using Scanning Tunneling Microscopy. ACS Chem Neurosci. 2011;2:281–287. [PMC free article] [PubMed]
31. Järvisalo J, Saris NE. Action of propranolol on mitochondrial functions--effects on energized ion fluxes in the presence of valinomycin. Biochem Pharmacol. 1975;24:1701–1705. [PubMed]
32. Freed AS, Garde S, Cramer SM. Molecular Simulations of Multimodal Ligand-Protein Binding: Elucidation of Binding Sites and Correlation with Experiments. J Phys Chem B. 2011;115:13320–13327. [PubMed]
33. Burra PV, Zhang Y, Godzik A, Boguslaw S. Global distribution of conformational states derived from redundant models in the PDB points to non-uniqueness of the protein structure. Proc Natl Acad Sci U S A. 2009;106:10505–10510. [PubMed]
34. Lill MA. Efficient Incorporation of Protein Flexibility and Dynamics into Molecular Docking Simulations. Biochemistry. 2011;50:6157–6169. [PMC free article] [PubMed]
35. Li DW, Showalter SA, Brüschweiler R. Entropy localization in proteins. J Phys Chem B. 2010;114:16036–16044. [PubMed]
36. Mobley DL, Dill KA. Binding of small-molecule ligands to proteins: “what you see” is not always “what you get” Structure. 2009;17:489–498. [PMC free article] [PubMed]
37. Butini S, Campiani G, Borriello M, Gemma S, Panico A, et al. Exploiting protein fluctuations at the active-site gorge of human cholinesterases: Further optimization of the design strategy to develop extremely potent inhibitors. J Med Chem. 2008;51:3154–3170. [PubMed]
38. Chiba S, Harano Y, Roth R, Kinoshita M, Sakurai M. Evaluation of protein-ligand binding free energy focused on its entropic components. J Comput Chem. 2012;33:550–560. [PubMed]
39. Jorgensen WL. The many roles of computation in drug discovery. Science. 2004;303:1813–1818. [PubMed]