Received date: December 06, 2016; Accepted date: January 05, 2017; Published date: January 08, 2017
Citation: Sorin EJ, Alvarado W, Cao S, Radcliffe A, La P, et al. (2017) Ensemble Molecular Dynamics of a Protein-Ligand Complex: Residual Inhibitor Entropy Enhances Drug Potency in Butyrylcholinesterase. Bioenergetics 6:145. doi: 10.4172/2167-7662.1000145
Copyright: © 2017 Sorin EJ, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Bioenergetics: Open Access
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 relaxation of these ensembles to conformational equilibria, binding modes for each inhibitor were identified. While classical models, which assume significant reduction in both 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.
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 [3-6]. These enzymes, which play physiochemically distinct roles both in neurotransmission and in cellular differentiation and development , 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 binding site (green), the acyl binding site (blue), and the peripheral anionic site (red).
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 [10-12]. 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 . The development of natural and synthetic cholinesterase inhibitors is thus a rapidly evolving field, and understanding both the physical interactions upon which proteinligand 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 , 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 approx. 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 [15-18] (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 , and which form the base units of intriguing dimeric inhibitors that show improved potency .
All-atom molecular dynamics simulations of native BChE  sans inhibitor, and of the protein in complex with each of these DAPP inhibitors, were performed using GROMACS 3.3  inside the Folding@Home distributed computing architecture . The ~60 kDa protein was modeled using the AMBER-03 force field of Duan, et al. . Inhibitors were modeled using the general AMBER force field (GAFF) with RESP charges derived at the 6-31G* level using RED Server . After initial docking of each inhibitor, sodium ions were randomly placed in a cubic 100 Å periodic box centered on the protein to establish electroneutrality. Solvation of this system with ~30,200 TIP3P  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  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 , 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-ofmass, (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 15-D vector was used as the basis for K-means clustering  to identify inhibitor binding modes.
|Project||RMSD* (Å)||Rg (Å)||Nhelix**||Nbeta**|
|Uninhibited||2.51 ± 0.25||23.02 ± 0.11||175.6 ± 6.9||86.3 ± 4.4|
|DAPP1 Ensemble||2.43 ± 0.24||23.02 ± 0.11||176.6 ± 6.9||84.9 ± 4.2|
|DAPP3 Ensemble||2.63 ± 0.30||23.06 ± 0.13||175.0 ± 7.0||87.5 ± 4.5|
|DAPP5 Ensemble||2.42 ± 0.21||23.02 ± 0.11||175.5 ± 6.6||86.2 ± 4.4|
Structural information provided in the final four rows of the table are averages over
all data in each ensemble following the 6.0 ns timepoint.
*RMSD values were calculated from the simulation starting structure of the enzyme after energy minimization, solvation, and solvent annealing.
**Number of residues in helical and beta conformations determined by DSSP , including β-sheet and β-bridge residues in the Nbeta category and all helical forms in the Nhelix category.
Table 1: Ensemble averaged structural metrics for BChE alone and in BChE-DAPP complexes.
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 each DAPP ensemble, representing moves from each binding mode i to each binding mode j after this equilibration period, were found to demonstrate both timeindependence 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 .
|KI* (mM)||1.7 ± 0.3||0.08 ± 0.02||0.006 ± 0.002|
|*Experimentally observed values for inhibition of BChE .|
Table 2: DAPP inhibitor binding properties.
Also shown in Table 2 are the experimentally observed inhibitor dissociation constants (KI) for these three inhibitors , which are the equilibrium constants for the undocking process BChEI → BChE+I, and which we can thus approximate based on the number of unbinding events observed in our ensemble simulations as KI=[BChE] [I]/[BChEI]. 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 inline with the tabulated experimental values 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. 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 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 residue compared to BChE sans inhibitor is shown for each BChE-DAPP complex, following the same color scheme used in the top panel. Dashed vertical lines in the background identify residues in the BChE binding pocket.
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 from these images for visual clarity, and residues in the binding pocket are colored as noted in the key. Binding modes are labeled from the bottom up, with mode 0 being the most populated (lowest binding free energy) and mode 23 being the least populated (highest binding free energy), and the binding mode free energy in kcal/mol (relative to the lowest energy mode) specified in the bottom right corner of each frame. All images are viewed down the ~20 Å deep BChE active site gorge from the same reference point and at approximately the same relative magnification.
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.  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 . 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 .
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 [33-35], 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 , 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 .
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 ligandreceptor 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].
Selected simulation movies are available at http://folding.cnsm. csulb.edu/BChE.php.
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.