**Cholinesterase Research**

Editors

**Ondrej Soukup Jan Korabecny**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Ondrej Soukup Biomedical Research Center, University Hospital Hradec Kralove Faculty of Military Health Sciences, University of Defense Hradec Kralove Czech Republic Jan Korabecny Biomedical Research Center, University Hospital Hradec Kralove Faculty of Military Health Sciences, University of Defense Hradec Kralove Czech Republic

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Biomolecules* (ISSN 2218-273X) (available at: www.mdpi.com/journal/biomolecules/special issues/ Cholinesterase Research).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-1798-8 (Hbk) ISBN 978-3-0365-1797-1 (PDF)**

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


Reprinted from: *Biomolecules* **2020**, *10*, 1556, doi:10.3390/biom10111556 . . . . . . . . . . . . . . . **151**

#### **Abdullah Al Mamun, Jana Maˇr´ıkov´a, Daniela Hulcov´a, Jiˇr´ı Janouˇsek, Marcela Safratov´a, ˇ Lucie Nov ´akov ´a, Tom ´aˇs Kuˇcera, Martina Hrabinov ´a, Jiˇr´ı Kuneˇs, Jan Kor´abeˇcny and Lucie ´ Cahl´ıkov ´a**

Amaryllidaceae Alkaloids of Belladine-Type from *Narcissus pseudonarcissus* cv. Carlton as New Selective Inhibitors of Butyrylcholinesterase

Reprinted from: *Biomolecules* **2020**, *10*, 800, doi:10.3390/biom10050800 . . . . . . . . . . . . . . . . **173**

#### **Gilbert Audira, Nguyen Thi Ngoc Anh, Bui Thi Ngoc Hieu, Nemi Malhotra, Petrus Siregar, Omar Villalobos, Oliver B. Villaflores, Tzong-Rong Ger, Jong-Chin Huang, Kelvin H.-C. Chen and Chung-Der Hsiao**

Evaluation of the Adverse Effects of Chronic Exposure to Donepezil (An Acetylcholinesterase Inhibitor) in Adult Zebrafish by Behavioral and Biochemical Assessments Reprinted from: *Biomolecules* **2020**, *10*, 1340, doi:10.3390/biom10091340 . . . . . . . . . . . . . . . **191**

**Faez Saleh Al-Hamed, Ola M. Maria, Jeff Phan, Ahmed Al Subaie, Qiman Gao, Alaa Mansour, Lina Abu Nada, Imane Boukhatem, Osama A. Elkashty, Simon D. Tran, Marie Lordkipanidz ´e, Zahi Badran and Faleh Tamimi**

Postoperative Administration of the Acetylcholinesterase Inhibitor, Donepezil, Interferes with Bone Healing and Implant Osseointegration in a Rat Model

Reprinted from: *Biomolecules* **2020**, *10*, 1318, doi:10.3390/biom10091318 . . . . . . . . . . . . . . . **217**

## **About the Editors**

#### **Ondrej Soukup**

Ondrej is the head of the Biomedical Research Center of the University Hospital Hradec Kralove and is focused mostly on the drug development process. He became the Associate Professor in Toxicology at the University of Defence in Hradec Kralove. He also obtained the PhD title in 2011 in this field. As a postdoctoral fellow, he was responsible for the development and testing of new therapeutics for Alzheimer's disease and organophosphorus poisoning. Besides the scientific activity, Ondrej is devoted to supervising students and teaching. The research focus can be described as a screening of biological properties of the new compounds, assessment of the blood-brain-barrier penetration, interaction with the cholinergic system and the determination of their cytotoxic effects; pharmacological and toxicological profile determination in vitro and in vivo, disinfectant and decontaminant agents and their efficacy.

#### **Jan Korabecny**

Jan is the deputy head at the Biomedical Research Center of the University Hospital Hradec Kralove and is mainly involved in the early stages of drug development. In 2020, he received an Associate Professor degree from Pharmaceutical Chemistry at the Faculty of Pharmacy in Hradec Kralove (Charles University). He also obtained his PhD title in 2012 in this field. His main area of interest is the development of novel therapeutics for Alzheimer's disease, cancer, tuberculosis, and organophosphorus poisoning. Besides the scientific activity, Jan is Editor-in-Chief of Military Medical Science Letters; he is supervising pregradual and PhD students. The research area is concentrated on drug development mainly applying different in silico techniques, and synthesis of small molecules as novel pharmacological entities. He is married and has two children.

## **Preface to "Cholinesterase Research"**

The scope of the Cholinesterase Research Special Issue was to provide a broad and updated overview of all the aspects that encompass cholinesterase research.

This collection of 10 scientific papers includes original as well as review articles focused on the cholinesterase structural aspects, drug design and development of novel cholinesterase ligands, but also contains papers focused on the natural compounds and their effect on the cholinergic system and unexplored effects of donepezil. The Special Issue is dedicated mostly to the scientific community in order to support research in the field of cholinesterases.

The editors would like to thank all authors for contribution to the special issue and for their hard everyday work in this field.

This work has been also institutionally supported by MH CZ –DRO (UHHK, 00179906) and by the Ministry of Defence of the Czech Republic "Long Term Development Plan –Medical Aspects of Weapons of Mass Destruction"of the Faculty of Military Health Sciences Hradec Kralove, University of Defence, Czech Republic.

> **Ondrej Soukup, Jan Korabecny** *Editors*

### *Editorial* **Cholinesterase Research**

**Jan Korabecny 1,2,\* and Ondrej Soukup 1,2,\***


Cholinesterases are fundamental players in the peripheral and central nervous systems. These serine hydrolases are presented by a two-membered family, namely acetylcholinesterase (AChE, E.C. 3.1.1.7) and butyrylcholinesterase (BChE, E.C. 3.1.1.8). Under physiological conditions, AChE terminates the action of acetylcholine at synapses. AChE is also implicated in the differentiation of embryonic stem cells, neuritogenesis, cell adhesion, synaptogenesis, activation of dopamine neurons, amyloid beta fiber assembly, haematopoiesis and thrombopoiesis, or regulation of glutamate-mediated hippocampal activity. Many compounds target to inhibit this enzyme in order to symptomatically counteract low cholinergic tone; however, irreversible AChE blockade may have fatal consequences. This phenomenon is typical for a class of highly toxic compounds—nerve agents and pesticides. The role of BChE is still extensively discussed; it plays an important role in cholinergic mediation, it contributes to neurogenesis, and has a detoxifying effect towards different xenobiotic drugs. It is also assumed that BChE overtakes the function of AChE in the case of malfunction or later stages of Alzheimer's disease (AD).

Based on the abovementioned, both AChE and BChE are considered as highly relevant targets in the field of medicinal chemistry. For neurodegenerative disorders such as AD, there is a strong consensus that AChE/BChE reversible inhibitors can, at least temporarily, alleviate the symptoms associated with the disorder, and enhance the cognitive performance of individuals. Other cholinesterase ligands, namely cholinesterase reactivators, typically endowed with strong nucleophilic function, can revert the irreversible action of organophosphorus compounds (nerve agents and pesticides). However, there are many other areas of research involving AChE and BChE, for example, pesticides; inflammation; and other neuronal disorders such as Lewy body dementia, Parkinson's disease, myasthenia gravis, and so on.

The scope of this Special Issue of Biomolecules, "Cholinesterase Research", was to provide a broad and updated overview of all the aspects that encompass cholinesterase research. The collection includes cholinesterase structural aspects, drug design and development, in vitro biochemical studies, animal studies, and computational approaches, all devoted primarily to cholinesterases.

Computational studies focus mostly on the structural and dynamical aspects of cholinesterases: first of all a comprehensive review of cholinesterase modelling and simulation was provided by De Boer and colleagues [1]. It analyses AChE/BChE structure and function using computer-based modelling and simulation techniques using different models of both enzymes. It also discusses key structural similarities in the active site gorges of the two enzymes, such as flexibility, binding site location, and function, as well as differences, such as gorge volume and binding site residue composition. Catalytic studies are also described, with an emphasis on the mechanism of acetylcholine hydrolysis by each enzyme and novel mutants that increase catalytic efficiency. The review also explores the inhibitory properties of several compounds currently approved by the FDA and other experimental drugs through Monte Carlo-based docking calculations and molecular dynamics simulations.

**Citation:** Korabecny, J.; Soukup, O. Cholinesterase Research. *Biomolecules* **2021**, *11*, 1121. https://doi.org/ 10.3390/biom11081121

Received: 27 July 2021 Accepted: 29 July 2021 Published: 30 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The study by Jo ´nczyk and colleagues [2] explores molecular mechanisms determining the efficiency and selectivity of individual oximes to reactivate AChE/BChE blocked by sarin and tabun. The study investigated the reactivation of AChE and BChE by selected oximes using molecular docking methods. It identified amino acids essential for effective reactivation and those responsible for the selectivity of individual oximes against inhibited AChE/BChE. The observation made herein can significantly contribute to support the search for new effective reactivators.

The paper of Lushchekina et al. [3] investigated the effect of different concentrations of sucrose on the protein and water dynamics in cholinesterases. It revealed a non-linear correlation with increasing sucrose concentration, i.e., first a decrease in the dynamics at 5 wt% followed by a gain at 10 wt% sucrose. The explanation of this phenomenon is that sucrose molecules interact with the surface of the protein and the entrance of the gorge at a lower concentration through the water layer, damping the motions at the surface, but increasing them inside the gorge. When increasing the sucrose concentration more, the sucrose molecules replace some of the water molecules at the surface, permitting again more water molecules to enter the gorge and opening simultaneously new pathways, among them the hypothesized backdoor to the gorge.

The study by Zueva and colleagues [4] is a kinetic study corroborated by molecular modelling simulation of human AChE inhibition by fluorinated acetophenone derivative, namely 1-(3-*tert*-butylphenyl)-2,2,2-trifluoroethanone (TFK). TFK was found to be a competitive type inhibitor reaching steady state inhibition slowly. It is speculated that after binding, TFK acylates the active serine, forming a hemiketal; the disruption of such complex, i.e., deacylation is slow. Modelling of interactions between TFK and AChE active site by QM/MM showed that the "isomerization" step of enzyme-inhibitor complex leads to a complex similar to substrate tetrahedral intermediate, a so-called "transition state analog", followed by a labile covalent intermediate. TFK can be classified as a slow-binding inhibitor with potential dual effect that could be of interest in palliative therapy of AD or protection of central AChE against organophosphorus compounds.

When it comes to a development of novel chemical compounds targeting cholinestarases both in the inhibitory and reactivation manner, the Special Issue contains the following experimental studies: the article by Konecny et al. [5] describes the design, synthesis and biological evaluation of a series of 15 novel fluoren-9-amine derivatives as dually active cholinesterase inhibitors and *N*-methyl-D-aspartate receptor (NMDAR) antagonists. The study builds on the concept of so-called multi-target directed ligands (MTDLs) that are believed to provide higher benefit compared to single-oriented drugs. The compounds under the study were initially in silico screened for CNS and oral availability, fitting all the prediction models used. Ongoing assessment of the biological profile included determination of the cholinesterase inhibition and NMDA receptor antagonism at the GluN1/GluN2A and GluN1/GluN2B subunits of NMDAR, along with a low cytotoxicity profile in the CHO-K1 cell line. Compounds were found to be highly selective BChE inhibitors with antagonistic activity on the NMDARs.

The group of Florian Nachon compared in vitro and in vivo efficacy, and toxicity of a hybrid tetrahydroacridine pyridinaldoxime reactivator, namely KM297, with pralidoxime [6]. The study revealed that blood–brain barrier crossing capacity of KM297 in vitro exceeds the permeability coefficient of pralidoxime twice. However, KM297 is also endowed with higher cytotoxicity, particularly on bone marrow-derived cells. Its strong cholinesterase inhibition potency seems to be correlated to its low protective efficacy in mice exposed to paraoxon. Ventilatory monitoring of KM297-treated mice by double-chamber plethysmography displayed toxic effects at the selected therapeutic dose.

Natural compounds and their effect on the cholinesterases are involved as well. Namely, the article by Amat-ur-Rasool and colleagues [7] screened methanolic extracts from seven commonly cultivated plants for their nutraceutical potential with particular emphasis on AChE/BChE inhibition and antioxidant capacity. The majority of extracts inhibited AChE and BChE, with henna and eucalyptus extracts highlighted as the most potent

ones. Moreover, all plant extracts were able to scavenge free radicals in a concentrationdependent manner, with eucalyptus being the most potent antioxidant.

The article by Al Mamun and colleagues [8] describes the isolation of thirteen known and three previously undescribed alkaloids of belladine structural type. Notably, significant human BChE inhibition was demonstrated by newly described alkaloids carltonine A and carltonine B, representing a new scaffold for generation of a novel structural type of BChE inhibitors with highly selective pattern over AChE.

Finally, the effect of clinically used donepezil has been investigated: the study of Audira and colleagues [9] is dedicated to donepezil, a currently approved drug for mildto-moderate stages of AD. The study exploits a zebrafish model to analyze potential adverse effects of donepezil on the short-term memory, behavioral and biochemical changes. Donepezil caused a slight improvement in the short-term memory of zebrafish and induced significant elevation in aggressiveness, while the novel tank and shoaling tests revealed anxiolytic-like behavior. The latter can be ascribed to alterations associated with an elevation of oxytocin and a reduction in cortisol levels in the brain. Thus, chronic waterborne exposure to donepezil can severely induce adverse effects on normal zebrafish in a dose-dependent manner.

In another article by Al-Hamed et al. [10], postoperative administration of donepezil on bone healing was studied in the group of Sprague–Dawley rats. After two weeks of donepezil administration, rats were euthanized, and their bones were analyzed by Micro-CT and histology, with the authors concluding that bone defects and implant osseointegration were significantly reduced compared to the saline-treated rats. Histomorphometric analysis pointed to lower immune cell infiltration in bone defects as the possible culprit for disrupted bone healing.

To conclude, this Special Issue describes important findings related to cholinesterases' physiological and pathological roles, their involvement in metabolic studies, and the influence of different modulators (inhibitors, reactivators) on their activity. All these findings may broaden the knowledge and the impact on clinical and pharmacological applications of different small molecules targeted to AChE/BChE.

**Funding:** The work has been supported also by the Long Term Development Plan—Medical Aspects of Weapons of Mass Destruction; of the Faculty of Military Health Sciences Hradec Kralove, University of Defence, Czech Republic and University Hospital Hradec Kralove (MH CZ-DRO) (UHHK, 00179906).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Review* **A Comprehensive Review of Cholinesterase Modeling and Simulation**

**Danna De Boer 1,†, Nguyet Nguyen 2,†, Jia Mao <sup>2</sup> , Jessica Moore <sup>3</sup> and Eric J. Sorin 1,\***


**Abstract:** The present article reviews published efforts to study acetylcholinesterase and butyrylcholinesterase structure and function using computer-based modeling and simulation techniques. Structures and models of both enzymes from various organisms, including rays, mice, and humans, are discussed to highlight key structural similarities in the active site gorges of the two enzymes, such as flexibility, binding site location, and function, as well as differences, such as gorge volume and binding site residue composition. Catalytic studies are also described, with an emphasis on the mechanism of acetylcholine hydrolysis by each enzyme and novel mutants that increase catalytic efficiency. The inhibitory activities of myriad compounds have been computationally assessed, primarily through Monte Carlo-based docking calculations and molecular dynamics simulations. Pharmaceutical compounds examined herein include FDA-approved therapeutics and their derivatives, as well as several other prescription drug derivatives. Cholinesterase interactions with both narcotics and organophosphate compounds are discussed, with the latter focusing primarily on molecular recognition studies of potential therapeutic value and on improving our understanding of the reactivation of cholinesterases that are bound to toxins. This review also explores the inhibitory properties of several other organic and biological moieties, as well as advancements in virtual screening methodologies with respect to these enzymes.

**Keywords:** acetylcholinesterase; butyrylcholinesterase; docking; molecular dynamics; hydrolysis; molecular recognition; catalysis; inhibition; reactivation

#### **1. Introduction**

The cholinesterase enzyme family has but two members: acetylcholinesterase (AChE) and butyrylcholinesterase (BChE). The former's primary biological purpose is regulating acetylcholine, a neurotransmitter, via hydrolysis at neuromuscular junctions, thus proving itself to be an essential component in the maintenance and performance of nervous systems. AChE, which has also been referred to as "true cholinesterase", is created in muscle, nerve, and hematopoietic cells and is considered to be one of the most efficient enzymes due to its rapid rate of catalysis [1]. AChE's plasma analog, BChE, previously referred to as "pseudocholinesterase", is produced in the liver and, unlike AChE, has been viewed as having a much more ambiguous biological purpose, as it was long believed to be vestigial [2]. BChE has a higher concentration in plasma than AChE, is present in many vertebrates, and tolerates several mutations, which has prompted the theory that BChE evolved from AChE to be a general detoxifier while still retaining some function in the process of neurotransmission [3]. This theory seems apt given the structural similarity of their binding sites, their shared affinities for certain substrates and ligands, and their sequence homology of approximately 65% [4]. Figure 1 depicts both AChE (bottom) and BChE (top) from a gorge-centric view and after a 90◦ rotation.

**Citation:** De Boer, D.; Nguyen, N.; Mao, J.; Moore, J.; Sorin, E.J. A Comprehensive Review of Cholinesterase Modeling and Simulation. *Biomolecules* **2021**, *11*, 580. https://doi.org/10.3390/biom11040580

Academic Editors: Jan Korabecny and Ondrej Soukup

Received: 21 March 2021 Accepted: 11 April 2021 Published: 15 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** Visualizations of (**a**) butyrylcholinesterase (BChE) (PDBID 1P0I) and (**b**) acetylcholinesterase (AChE) (PDBID 1B41) in white ribbon mode with residues in notable binding sites shown as semitransparent van der Waals surfaces colored according to the key. Top panels present views directly into the active site gorges with structures in the center panel rotated 90◦ about the vertical axis. Bottom panels present magnified views of the binding site regions of (**c**) BChE and (**d**) AChE with key residues labeled.

One difference between the two enzymes is the relative size of the binding pocket, with BChE and AChE having approximate gorge volumes of 1500 Å<sup>3</sup> and 1300 Å<sup>3</sup> , respectively [5]. X-ray crystallography of *Torpedo californica* (*Tc,* pacific ray) AChE revealed the enzyme to have a deep hydrophobic gorge with residues that stabilize substrates in the pocket [6], as well as a bottleneck region in the active site [7] that narrows to approximately 4 Å in width [8]. Common models of AChE, including human, mouse, and *Torpedo californica* (*Tc*AChE), demonstrate conserved active sites, save for a few residues that participate in ligand binding [9], and both have negative surface potentials that become more negative deeper within the gorge. This negative potential, which is high near the catalytic site at the "bottom" of the gorge, seems to have evolved to facilitate electrostatic attraction of positively charged choline substrates by both enzymes [8].

Although BChE is structurally similar to AChE, 6 of the 14 aromatic amino acids that line the active site gorge in AChE are substituted with aliphatic residues in BChE [10]. In particular, the substitution of Phe288 and Phe290 in *Tc*AChE with the smaller Leu286 and Val288 in human BChE lead to conformational changes that result in a deeper gorge in BChE, thereby allowing BChE to interact with, and potentially hydrolize, a much wider range of substrates and inhibitors than AChE [11]. BChE is thus characterized as the promiscuous, or non-specific, bigger sibling to the smaller and much more specific AChE.

The tremendous growth and improvement in computational resources and modeling techniques over the past few decades have led to an exponential increase in computational studies of biomolecular systems. In 2003, the Protein Data Bank (PDB) launched its online presence, thereby making myriad biomolecules and complexes available for structural and computational research, and by 2016 there were 178 AChE structures available to the public, many of those including bound substrates or inhibitors of medical and pharmacological significance [9]. In fact, some of the earliest and most significant contributions of cholinesterase models were reported by Sussman et al., who performed an X-ray analysis of AChE at 2.8 resolution in 1991 [12], and Nicolet et al., who examined several crystal structures of BChE in 2003 [13]. Modeling and simulation-based studies of the cholinesterase enzymes have advanced in tandem with this growth.

The application of molecular dynamics (MD) simulations, Monte Carlo (MC) based docking calculations, and more sophisticated quantum mechanical/molecular mechanics (QM/MM) simulations have proven highly insightful in probing cholinesterase structure and activity. Notable areas of interest include the mechanisms of cholinesterase catalysis, reversible and irreversible inhibition of both enzymes to manage Alzheimer's Disease (AD) and other ailments, BChE-specific inhibition, and the reactivation of phosphorylated cholinesterases, a process that normally follows nerve agent attacks or pesticide poisoning [14]. Figure 2a,b depict the common paths of substrate or ligand binding and phosphorylation, respectively.

Most therapeutic treatments of AD, and other illnesses to which the cholinesterases have been linked, are reversible inhibitors that form non-covalent molecular recognition (MR) complexes with cholinesterases and can leave the active site, thus existing in an equilibrium between bound and unbound states characterized by K<sup>I</sup> and/or IC<sup>50</sup> values. The three mechanisms that reversible inhibition can follow are competitive, noncompetitive, and uncompetitive, as depicted in Figure 2a. During competitive inhibition, the substrate and inhibitor are competing for the same binding site, and the binding of a competitive inhibitor within the active site blocks entrance of the substrate, thereby hindering formation of the enzyme substrate complex. By contrast, uncompetitive inhibitors function by binding to the enzyme-substrate complex to prevent product formation, and noncompetitive inhibitors can bind to either the enzyme or the enzyme-substrate complex to regulate catalytic activity.

Many organic compounds, including organophosphates (OPs), are categorized as irreversible inhibitors, which covalently bond to residues in the gorge and thus cannot leave the active site. In the cholinesterases, OPs covalently bind to the serine residue of the catalytic triad. Irreversibly inhibited cholinesterases, however, can be reactivated via various pathways with compounds such as oximes, as illustrated in Figure 2b. Oximes reactivate an inhibited cholinesterase to its native structure by nucleophilic substitution of the phosphorylated serine. If untreated, the inhibited enzyme can undergo the dealkylation process, also known as aging, where the loss of the second leaving group produces an oxyanion on the phosphoryl group.

**Figure 2.** (**a**) Depiction of the enzyme catalytic mechanism (black) and reversible inhibition mechanisms including competitive inhibition (blue), uncompetitive inhibition (red), and noncompetitive inhibition (both blue and red). (**b**) Representation of irreversible (covalent) inhibition by organophosphorus ligands and the reactivation mechanisms to revert the aged enzyme (blue) back to the inhibited complex (green). L and R represent the first and second leaving groups of the organophosphate, respectively, following the paths of inhibition and reactivation.

Aged cholinesterase is highly stable due to the strong electrostatic interactions between the oxyanion and the positively charged catalytic histidine. Known reactivators of OP-inhibited AChE, such as oximes, are ineffectual against aged AChE [15]. Though a successful attempt to reactivate aged AChE by a class of compounds called "quinone methide precursors" (QMPs) was reported by Hadad and coworkers [16], the reaction between QMPs and aged AChE was rendered too slow to be useful [17]. The present review will examine many studies generally focusing on these areas of interest and, unless otherwise specified, all amino acid numbering used below will refer to human cholinesterase models.

#### **2. Structure and Dynamics**

While understanding the structure, function, and behavior of cholinesterase binding sites seems a practical starting point for modeling studies, early researchers initially wanted to understand how substrates and inhibitors entered the active sites of these enzymes. Molecular dynamics (MD) simulations of human BChE lasting for 5 and 10 ns indicated that inhibitors can access the binding subsites in the catalytic cavity due to the highly flexible entrance (or "mouth") of the gorge, a portion of which is formed by the flexible omega loop (Ω-loop) region, as well as the peripheral aromatic site (PAS, formerly mistakenly called the peripheral anionic site). In these simulations, the Asp70 residue in PAS showed significant deviation from the crystal structure with root-mean-square deviation (RMSD) values of 2 to 6 Å [18]. AChE has also been shown to experience such fluctuations at the gorge entrance, with a similarly flexible Ω-loop region that is thought to increase enzyme specificity by making it more difficult for large molecules to enter the AChE gorge without hindering the productivity of the enzyme [19].

Apart from the fluctuations at the mouth of the gorge, AChE also experiences what is known as bottleneck fluctuations, with the bottleneck being the narrowest part of the

gorge. These fluctuations, which have come to be known as the "breathing" of the enzyme, can help inhibitors or substrates move from the surface to deeper regions of the gorge [20]. Cheng et al. recently defined this breathing by monitoring the varying distance between the Cε<sup>2</sup> atom of Phe330 (CBS) and the O<sup>H</sup> atom of Tyr121 (PAS) in TcAChE, which suggested that a number of subdomains within the enzyme, particularly the Ω-loop, contribute to modulating the size of the gorge bottleneck [21]. A comprehensive comparison between 47 crystal structures of AChE (in its native form and in complex with small molecules), as well as a 20 ns simulation of *Tc*AChE, provided by Sussman and coworkers, suggested that the 14 aromatic residues lining the AChE gorge and creating over half of the gorge surface area contribute greatly to the overall flexibility of the enzyme, the observed bottleneck breathing motions, and the resulting ability to perform its catalytic function [22].

Some of these aromatic residues play significant roles in primary subsites within the cholinesterase gorges, such as the catalytic active site (CAS) and the peripheral aromatic site (PAS) [22]. While these sites, which are addressed in more detail below, are integral to cholinergic activity, cholinesterases do not only perform cholinergic functions. As discussed by Chinnadurai et al., the aryl acylamidase activity (AAA) of AChE, which also involves hydrolysis, only requires the CAS and does not interact with the PAS at all [23]. This has prompted the theory that AAA substrates enter from a side-door into the enzyme, rather than via the mouth of the gorge and, indeed, researchers have suggested that there are a number of doors through which substrates can enter the gorge including a back door [24], an Ω-loop door, and the suggested side door [23]. It was suggested that the side door may open more frequently than the other doors to mediate AAA activity, and MD simulations of side door probing emphasize the importance of hydrophobic interactions, hydrogen bonding, and water mediated interactions ("water bridges") in moving the substrate towards the CAS [23]. To be sure, simulations of AChE in explicit solvent sans substrate [25], as well as analyses of *Tc*AChE crystal structures in its native and several inhibited forms [26], have underscored the importance of the presence of molecular water in enzyme structure.

At higher concentrations, dimerization and the further dimerization of dimers to form tetramers is known to affect the structure and function of cholinesterase enzymes [27]. While MD simulations suggested that two of the four binding sites in tetramerized cholinesterases are sterically blocked, thereby becoming less active [28], as reflected by a 15% decrease in catalytic activity [29], a recently elucidated CryoEM structure of hBChE shows distinct structural variance from the simulated tetramer, with the active site gorges being fully solvent accessible [27]. Tetramerization, however, increases the half-life of the enzymes, which is a desirable result when cholinesterases, and BChE in particular, are being used to counteract drug overdoses. For example, the addition of proline-rich attachment domains (PRADs) to BChE increases tetramer stability, leading to an extended circulation time [30]. Interestingly, glycosylated models of BChE increase the enzyme's flexibility and half-life without hindering its ability to bind to glycans, which cannot be said of all therapeutic protein targets [31]. On the other hand, cholinesterase phosphonylation, or the irreversible binding of organophosphate to the active site, which will be discussed in more detail in the organophosphate inhibition section below, severely restricts the flexibility of both AChE and BChE, as reported by Bennion et al. [32]. Experimental findings of AChE covalently bound to the nerve agent soman agree that OP-poisoned AChE is significantly stiffer [33].

#### *2.1. Important Binding Sites*

There has been substantial past effort to study the sites responsible for molecular recognition (MR) and binding affinity within the active site gorge of both enzymes, and it is clear that those binding sites, specific chemical subsites within each active site gorge, mirror each other and perform similar functions respective to each enzyme. While it is important to note that these binding sites have been examined experimentally with X-ray and kinetic studies, including a recent study by Rosenberry et al. [34], the present review

focuses on the unique perspective provided by computational investigations. As expected, one of the most important sites for both cholinesterases is the catalytic active site (CAS), and the peripheral aromatic site (PAS) also plays an indispensable role in cholinesterase or ligand binding, while the Ω-loop (OML), acyl binding site (ABS), and oxyanion hole (OAH) sites are more essential in contributing to binding affinity and complex stability. Alvarado et al. have provided a method of succinct graphical tabulation of BChE-ligand contacts and interactions, referred to as contact tables, that include these five sites and additional protein residues of interest [35], as discussed below. A detailed analysis of these binding sites is provided here in the same order that they are encountered by substrates and inhibitors upon entering and moving into the gorge.

#### 2.1.1. Peripheral Aromatic Site

The peripheral aromatic site (PAS, red in Figure 1) is located near the mouth of gorge [36] and plays a prominent role in substrate and ligand binding [37]. For decades, peer-reviewed studies have used PAS to denote the "peripheral anionic site". In recent years, however, the aromatic properties of this binding site that are vital to cholinesterase function have driven the community to instead refer to this region as the "peripheral aromatic site". Important amino acids in the PAS of AChE include serine, tyrosine, aspartic acid, and tryptophan [36,38], while notable PAS residues in BChE include asparagine, aspartic acid, glutamine, serine, and tyrosine [35], highlighting the polar, negatively charged, and electron-rich nature of residues in this site. As previously mentioned, one distinction between the cholinesterases is the aromatic nature of the residues surrounding the PAS in AChE, which is more aliphatic in BChE [10]. The PAS makes contact with many loops and secondary structural elements at the surface of the protein, including the Ω-loop, which contributes to the needed flexibility discussed above. Although steric and electrostatic interactions may slow the catalytic efficiency of AChE, the PAS is valuable for trafficking ligands into the gorge [39], particularly positively charged species such as cholines. MD simulations have emphasized the importance of cation-π interactions, which stabilize the ligand at the rim of the gorge entrance prior to entering the gorge [40], and it has been proposed that non-cholinergic activity of the PAS could include the deposition of amyloids, adhesion to cells, and outgrowth of neurites [41].

#### 2.1.2. Acyl and Choline Binding Sites

Once a ligand has entered the gorge, the acyl and choline binding sites (ABS and CBS, shown as blue and green in Figure 1, respectively), which are located near the catalytic triad, assist in positioning the ligand for catalysis. The ABS and CBS are hydrophobic regions composed primarily of tryptophan, tyrosine, and phenylalanine in human AChE. Tyr337 in the choline binding site of AChE is replaced by Ala328 in that of BChE; Phe295 and Phe297 in the acyl binding site of AChE are replaced by Leu286 and Val288, respectively, in BChE. The replacement of aromatic residues in the ABS and CBS of BChE enable it to bind larger substrates and inhibitors than AChE [42]. In addition, the ABS and CBS are largely responsible for the specificity of these enzymes and are thus primary targets studied when synthesizing inhibitors such as imidazole or pyridine derivatives [43].

#### 2.1.3. Catalytic Active Site

The catalytic active site (CAS, yellow in Figure 1) is located approximately 20 Å deep at the bottom of both the AChE and BChE gorges [12,24] and is made up of serine, glutamic acid, and histidine residues, prompting the name "catalytic triad" [24,35]. The CAS is surrounded by numerous aromatic and acidic residues [44] and is observed to engage in shorter hydrogen bonds in crystal and NMR structures than observed in simulation [45,46]. More importantly, the CAS is highly conserved [44], emphasizing historically vital biological roles of these enzymes and their cholinergic activities. QM/MM simulations at the MP2(6-31 + G\*) level reveal a potential energy barrier of 10.5 kcal/mol, which agrees with experimental data [47], and MD simulations of AChE bound to acetylcholine (ACh)

show that ACh stabilizes the CAS and improves the binding ability of the peripheral aromatic site [48].

It has been postulated that a back-door exists in AChE, just behind the CAS and controlled by a tryptophan residue, which was theorized after a single water molecule exited the active site gorge from a direction opposite that of the gorge entrance in an MD simulation [49]. This "back door" was later thought to open three to four Å wide such that catalysis products could exit the gorge of the enzyme without blocking the gorge entrance, and thus contributing to a high catalytic rate [50]. Aromatic residues surrounding the CAS histidine also largely influence the productivity and efficiency of the enzyme [51], which decreased approximately 600-fold when disrupted or replaced by aliphatic side chains [52]. More recently, Xu et al. used MD to study *Tc*AChE and, from 27 of their 40 trajectories, observed thiocholine to frequently exit the active site via a back-door created by cooperative motions of CAS residue Trp86 along with Val132 and Gly448 [53], for which previous experimental support was noted [54,55].

Indeed, other mutations in or near the CAS are known to have effects on the structure, and subsequent function, of the enzyme. This research has naturally focused on, and is more applicable to, BChE due to its much greater natural affinity for mutations than AChE [3]. For example, prolonged use of muscle relaxers led to the discovery of the "silent phenotype" in which an alanine is mutated to a valine near the CAS of BChE. This mutation was studied in silico and observed to severely disrupt interactions between the histidine and serine in the catalytic triad [56], leading to a dysfunctional CAS, regardless of the inhibitor, and increases in the volume of the enzyme, indicating that this may be a pre-denaturation state [57].

The mutation of a nearby alanine in silico to cystine in wild-type (WT) BChE causes the histidine in the catalytic triad to flip, an event that is largely guided by local water molecules [58]. This man-made mutation, while possibly slowing the speed of binding, ultimately still allows for substrate binding to the active site; the naturally occurring mutation of that alanine to aspartic acid, however, is claimed to be catalytically inactive due to strong disruptive interactions between aspartic acid and the CAS histidine residue [59]. Both mutations showcase the possible hysteretic behavior of BChE, or its reliance on past-states, which can likely be attributed to its toxicological and pharmaceutical functions [58,59]. For instance, a man-made mutant of BChE was recently modeled and examined by Masson and coworkers, using QM/MM and Markov state analysis, and was determined to be a template for future investigations into organophosphate hydrolase functions [60].

#### 2.1.4. Oxyanion Hole

The oxyanion hole (OAH, orange in Figure 1), made of two glycines and one alanine [35,47], is generally a two-pronged site in many proteases and hydrolases; in the case of AChE, and subsequently BChE, the OAH is a three-pronged hole [61]. Early MD simulations of AChE phosphonylation suggested that the OAH exerts a pulling force on leaving groups during the alkylation step [45], and the OAH is known to lower the energy barrier for ACh hydrolysis in both cholinesterases [62]. QM/MM simulations exhibited consistent, tight hydrogen bonding between the OAH and the carbonyl carbon of the substrate, suggesting that the OAH facilitates stabilizing interactions in intermediate and transition states [61].

#### 2.1.5. Ω-loop (Omega Loop)

The Ω-loop (OML, purple in Figure 1), consisting of a series of nearly 30 residues, is located along one side of the active site gorge wall. In the presence of a substrate or inhibitor, the Ω-loop plays an important role in modulating enzyme "breathing", regulating the size of the gorge, and thereby enabling the passage of the ligand to the active site [21]. The OML undergoes conformational changes, such as gorge enlargement, through torsional motion and segmental fluctuations [63,64]. Unregulated motions and decreased electrostatic interactions, however, can significantly decrease the binding affinity of these

enzymes, such as the case of the atypical mutation from Asp70 to Gly70 in the OML of BChE, which Masson et al. reported could increase the *K*<sup>m</sup> values 10- to 100-fold [65]. Moreover, the Ω-loop is speculated to facilitate an alternative entrance, the proposed sidedoor model noted above, to the active site of the enzyme; MD simulations performed by Wiesner et al. suggested that protonation of the AChE active site leads to conformational changes within the Ω-loop at Asn87 and Glu84 that result in the opening of this side door [66]. A similar observation was recently reported by the Rydzewski laboratory for *Tc*AChE, where opening and closing of this side door due to the displacement of the OML favored alternative dissociate routes of the substrate and inhibitor [67]. To put these computational results into perspective, a number of experimental investigations into the omega loop and backdoor of AChE from various species have suggested that the back door opens in some cases [50,68,69], but also that this opening is likely not relevant functionally [70].

#### **3. Catalysis**

The primary catalytic function performed by AChE is hydrolysis of the neurotransmitter acetylcholine (ACh) [1], and molecular modeling and simulation studies have provided insight into this process that cannot be easily gleaned from experimental efforts. For example, MD simulations using an ab initio QM/MM potential were conducted to map the reaction mechanism of AChE with ACh, pointing to a mechanism with two main processes: acylation and deacylation. In the first step of acylation, the system must overcome the initial free energy barrier of 12.4 kcal/mol, during which the oxygen atom of Ser203 performs a nucleophilic attack at the carbonyl carbon of acetylcholine, with a synchronous proton transfer from Ser203 to His447, resulting in the first tetrahedral intermediate [71]. This intermediate is stabilized by local hydrogen bonds from the oxyanion region and electrostatic interactions with the glutamic acid in the catalytic triad [47]. The second step in acylation has an additional free energy barrier of 1.9 kcal/mol and is characterized by the proton transferring to the leaving group of ACh, resulting in bond cleavage and choline release from the pocket [71].

The next step in this process, deacylation, also occurs in two steps, the first of which includes interaction between a water molecule and the carbonyl carbon of the acetyl group in acetylcholine, leading to a second tetrahedral intermediate that is stabilized by the Gly121, Gly122, and Ala204 residues of the oxyanion hole. The proton then transfers to the acetylserine oxygen atom from His447, yielding the products acetic acid and AChE, with the initial free energy barrier of deacylation higher than that of acylation at 17.5 kcal/mol and thus predicted to be the rate-limiting step [71].

Modeling efforts by Chen et al. have shown BChE to have a similar reaction pathway for hydrolysis of ACh and acetylthiocholine, with two-step acylation and deacylation processes [72,73]. Like AChE, residues in the oxyanion hole of BChE are integral for stabilization of the intermediate and catalytic function [73]. However, the acylation and deacylation processes for BChE were predicted to have free energy barriers of 13.8 and 11.9 kcal/mol, respectively, indicating that the rate-limiting step for BChE is not the deacylation step, as predicted for AChE, but rather the acylation step [72]. Figure 3 displays free energy profiles for the hydrolysis of ACh by both AChE and BChE to allow for a side-by-side comparison [71,72] of these results.

Given the promiscuous nature of BChE and its circulation in plasma, this enzyme has a much broader natural variety of substrates than AChE [2,3]. For example, BChE is one of the primary enzymes to hydrolyze heroin and produce its most active form, 6-monoacetylmorphine. Computational efforts by Zhan and coworkers showed this process to follow the two-step acylation and deacylation scheme outlined above [74]. In contrast, work by the same group to model BChE hydrolysis of ghrelin, the hunger hormone, showed a single-step acylation process [75]. Interestingly, Suarez et al. modeled BChE hydrolysis of butyrylcholine, a synthetic molecule that mimics acetylcholine and for which BChE is named [3], and found that the presence of glycerol or another butyrylcholine

in the active site pocket stabilizes the intermediate product after transition state 2 in the deacylation step, which is a complex of BChE with butyric acid [76].

**Figure 3.** (**Top**) Free energy profiles for hydrolysis of ACh by AChE including the (**a**) acylation and (**b**) deacylation steps [71]. (**Bottom**) Free energy profiles for hydrolysis of ACh by BChE showing analogous (**c**) acylation and (**d**) deacylation steps [72].

Cocaine is another substrate of interest for BChE catalysis studies, as increasing the catalytic efficiency of BChE for cocaine hydrolysis can be an effective method to treat overdoses, making BChE mutants and transition states important focal points for study. Rate-determining steps can be made faster with residue mutations [77] by changing local interactions and increasing substrate stability [78]. For example, computational mutations of non-active residues in the BChE active site gorge were predicted to increase catalytical efficiency as much as 1000-fold by strengthening hydrogen bonds [61]. The Ala199Ser/Ser287Gly/Ala328Trp/Tyr332Gly BChE mutant achieves this by increasing the strength of hydrogen bonds in the first transition state of BChE-cocaine catalysis and lowering the energy barrier [79].

Free energy perturbation (FEP) simulations have allowed researchers to see the deviation in free energy barriers at transition states for different BChE mutants and led to the discovery of a mutant, Ala328Trp/Ala199Ser/Phe227Ala/Glu441Asp/Ser287Gly, that is around 1800-fold more efficient than wild-type BChE [80]. Further work by Zhan and coworkers using FEP simulations showed the Ala328Trp/Tyr332Gly/Ala199Ser BChE mutant to have the potential to greatly increase catalytic efficiency and to thus serve as a potential means for exogenous therapy [81]. Using QM/MM simulations, Zhan's team revealed two transition states for BChE-cocaine binding involving deformation of the nonprereactive complex and formation of the prereactive complex [77,82]. Indeed, efforts by the Zhan laboratory have resulted in a prodigious quantity and breadth of computational studies of cholinesterase catalysis in the 21st century.

#### **4. Inhibition**

We now turn to the area around which the majority of cholinesterase research has been dedicated over the past two decades: inhibition of one or both enzymes by various compounds, presented below in distinct chemical groupings including pharmaceuticals, narcotics, organophosphates, other organic species, and biological agents and salts. FDAapproved pharmaceutical inhibitors and their derivatives generally serve to treat the symptoms of AD. Due to the difference in size of the two enzymes, many drugs experience a wide range of selectivity and target multiple binding sites. As previously mentioned, the biological roles that BChE plays are much more ambiguous than those of AChE. However, hydrolysis of narcotics, such as nicotine and cocaine, has been identified as a possible role of BChE and, as noted above, this enzyme has been a target to treat narcotics overdoses.

Organophosphate inhibitors have proven to have a wide field of study with myriad applications, and many phosphate-based molecules serve as potential cholinesterase inhibitors for disease management, while also being known for their highly toxic roles as irreversible inhibitors in nerve agents and pesticides. Finally, there are myriad compounds being researched, many aimed specifically at AD and other human ailments, that are still in early phases of research and not yet FDA-approved. Due to the vast number of unique inhibitors studied over the past two decades, many have been classified below under the broad category of other organic species, which are further organized according to prominent functional groups. Prior to this current effort, Anand and Singh reviewed different classes of cholinesterase inhibitors including tacrine, donepezil, rivastigmine, galantamine, xanthostigmine, para-amino-benzoic acid, coumarin, flavonoid, and pyrrolo-isoxazole analogues in their 2013 article [83]. Figure 4 presents molecular structures for many of the pharmaceuticals, narcotics, nerve agents, and related inhibitors mentioned above and detailed in the sections below.

#### *4.1. Pharmaceuticals*

#### 4.1.1. Tacrine and Derivatives

Tacrine, or 1,2,3,4-tetrahydroacridin-9-amine, was a commonly-used drug to treat AD (under the brand name Cognex) that was the first FDA-approved cholinesterase inhibitor, but was later discontinued due to liver toxicity [84]. Given this history, it is used as a comparison for other drug studies [85], as well as the subject of study for potential chemical derivatives. Intensive 3D-QSAR (quantitative structure-activity relationship) studies, molecular docking analyses, and MD simulations of 60 tacrine-based inhibitors bound to *Tc*AChE have identified key residues involved in binding to be Tyr70 (PAS), Trp84 (CBS), Tyr121 (PAS), Trp279 (PAS), and Phe330 (CBS) across various binding sites [86].

Tacrine derivatives have been investigated using docking calculations. In a survey of racemic tacrine derivatives in complex with AChE, Maalej et al. found the derivative 4-(13-amino-10,11,12,14-tetrahydro-9H-benzo[5,6]chromeno[2,3-b]quinolin-14-yl)phenol to be four times more effective in inhibitory activity than tacrine [87]. Tacrine-carbazole hybrids in complex with both AChE and BChE were also investigated via docking by Thiratmatrakul et al., who not only found these hybrids to exhibit a preference for BChE, but also showed them to have a potential for ABTS radical scavenging activity [88]. In the early 2000s, in situ click-chemistry led to the discovery of potent, femtomolar range tacrine derivatives [89], which have since been studied in complex with mouse, *Drosophila melanogaster*, and *Tc*AChE through MD simulations [90]. More recently, the Richardson laboratory employed docking and quantum characterization, which revealed that tacrine derivatives with spacers, such as pentylaminopropene and pentylaminopropane, can increase inhibitor specificity to target BChE [91].

As expected, MD simulations can provide more detailed insight into tacrinecholinesterase complexes, especially when coupled with experimental observations. For example, MD simulations allowed Decker and coworkers to explore interactions between ring-opened and ring-closed cyclohexen-like rings of tacrine derivatives [92], as well as to provide IC<sup>50</sup> values for indole-3-acetic acid (IAA)-tacrine dual AChE/BChE inhibitors [93]. Other studies using only computational methodologies have focused on tacrine-cholinesterase complex dynamics. While recent QM/MM simulations of tacrine bound to AChE reported by Nascimento et al. demonstrated that van der Waals forces play

just as important a role as electrostatics in binding and stabilization of the complex [94], MD simulations of BChE in complex with tacrine reported by Wan et al. emphasize the importance of protonating Glu197 near the catalytic triad, which stabilizes participation of a localized water molecule and leads to preservation of the His438 residue [95].

**Figure 4.** Chemical structures of various narcotic, nerve agent, and pharmaceutical cholinesterase inhibitors discussed in this review.

Interestingly, MD simulations of a survey of tacrine-huprine heterodimers revealed that derivatives show potential for AD and prion disease treatment because they both inhibit the PAS and CAS in AChE while reducing β-amyloid and prion peptide aggregation [96]. Other dimers that exhibit dual-site inhibitory activity with the cholinesterase enzymes are the bis(7)tacrine derivative with S-allylcysteine and cystamine, both of which can also serve as antioxidants [97,98]. The hybrid of tacrine and quercetin, dubbed Tac-Quer in a 2019 study by Habibpour et al., allows tacrine to act as an inhibitor to either AChE or BChE, while quercetin, a well-known metal scavenger, seeks Zn2+, Cu2+, and Fe2+ cations in trace quantities in brain plaque [99].

#### 4.1.2. Galantamine and Derivatives

Galantamine, or (4aS,6R,8aS)-5,6,9,10,11,12-hexahydro-3-methoxy-11-methyl-4aH- [1]benzofuro[3a,3,2-ef][2]benzazepin-6-ol, is another pharmaceutical that targets AD, sold under the brand name Razadyne in the United States [84]. It can be isolated from the *Amaryllidacease* family of plants [100] and, like tacrine, is used as a reference compound in drug discovery [4]. Examination of the X-ray structure of the *Tc*AChE-galantamine complex hinted that galantamine interacts with residue Trp84 at the choline binding site and residues Phe288 and Phe290 at the acyl binding site [38]. Recent docking calculations, however, suggest that galantamine inhibits AChE at the base of the active site gorge [101], and MD simulations point to the importance of hydrogen bonds with water for inhibitor stability in the gorge [102]. While galantamine is larger than the AChE gorge entrance observed in crystal structures, a recent experimental-computational collaboration by Roca et al. illustrated that it can enter the active site pocket of the enzyme following reorientation of the PAS to traffic the ligand inside, as mentioned previously [39].

#### 4.1.3. Donepezil and Derivatives

Donepezil, or 1-benzyl-4-[(5,6-dimethoxy-1-indanon-2-yl)methyl]piperidine, sold under the brand name Aricept [84], also serves as a reference compound in drug discovery studies [102], and was shown in a recent MD study by the Treptow laboratory to act as a mixed competitive and non-competitive inhibitor that interacts strongly with the PAS, ABS, and CAS regions of AChE [103].

QSAR examinations into cholinesterase inhibition have produced handfuls of donepezil derivatives [104,105], indicating that the parent compound is structurally favored to be an inhibitor. Docking simulations of pyridonepezil and quinolinodonepezil derivatives were performed with both cholinesterases, with quinolinodonepezil derivatives proving much less effective at inhibiting human AChE than pyridonepezil derivatives [106]. Marco–Contelles and coworkers took a multi-pronged approach to study donepezil-pyridyl hybrids, which proved to inhibit both cholinesterases at the PAS and CAS, and suggested that the N-alkyl bridge could be used to selectively enhance AChE inhibition by such pyridine-based derivatives [104]. Docking studies by Al-Rashid and Hsung further suggest that the E-ring in the donepezil-like (+)-arisugacin A compound can play a crucial role in binding to AChE [107], and Rahman et al. employed DFT and docking to demonstrate that halogenated derivatives of donepezil, including fluorine and chlorine groups, also show AChE inhibiting potential [108].

In a computational-experimental collaborative study focusing on a series of N-substituted amine derivatives [84], docking and MD simulations were performed using the AChE crystal structure from the pre-formed AChE-donepezil complex and some compounds were observed to mimic the binding pose (position and orientation) of donepezil [109]. In another more recent study, well-tempered metadynamics (WTMtD) simulations of AChE in complex with donepezil by Ghosh et al. showed the protein-ligand complex to increase the ordering of water molecules around Ser203 of the CAS, which discourages ACh from interacting with the active site [110]. Hybrids of donepezil's benzylpiperidine moiety connected via an oligomethylene linker to an indolyl propargylamino moiety were examined via MD simulations and identified as dual-site binding cholinesterase inhibitors [111]. Docking and

MD simulations by Yekta et al. suggest that glycated-AChE (a glycine and AChE complex) poses a challenge for donepezil binding due to the rearrangement of Trp286 and Tyr341 that block this inhibitor from entering the binding cavity [112].

#### 4.1.4. Rivastigmine and Derivatives

Exelon is the brand name for rivastigmine, or [3-[(1*S*)-1-(dimethylamino) ethyl]phenyl] *N*-ethyl-N-methylcarbamate [84]. Docking and MD simulations conducted recently by Ali et al. show that, due to the presence of rivastigmine, TcAChE undergoes carbamylation [102]. Rivastigmine and numerous conformationally restricted analogs were studied by Bolognesi et al. using Monte Carlo calculations, which suggested that the carbamic N-alkyl chain has a more negative effect on binding to AChE than BChE due to the larger acyl binding site present in BChE [113]. Another much more recent study by Wang et al. focusing on a series of chalcone-rivastigmine hybrids support this observation, with results from MD simulations showing rivastigmine hybrids to bind to BChE more easily [114].

#### 4.1.5. Quinazoline and Derivatives

Quinazoline is another pharmaceutical that has been considered as a potential AD treatment. In fact, a number of derivatives have been FDA-approved as anti-cancer and anti-tumor drugs, including Gefitinib, Erlotinib, Vandetanib, Lapatinib, and Afatinib [115]. Although, of the two cholinesterases, quinazoline derivatives seem to bind more effectively to AChE [116], the addition of alicyclic groups to quinazoline analogs increases the binding affinity towards BChE, as these groups bind more effectively to the PAS [117]. Homobivalent quinazolinimes are also derivatives that bind more closely with BChE; here, in docking simulations, the homobivalent quinazolinimes engage with BChE with π-interactions that are absent in the AChE complex [118]. However, in a study by the Decker laboratory tricyclic and tetracyclic quinazoline derivatives led to an "inverted binding mode" with the aliphatic amine in the center [119] and Daoud et al., who used a multi-pronged computational approach to study cholinesterase inhibition, found that pyrazinamide derivatives exhibit strong hydrogen bonding with Tyr121 in *Tc*AChE and Tyr332 in BChE, suggesting that these derivatives are ChE effective inhibitors [120].

#### 4.1.6. Coumarin

Coumarin, or 2H-1-benzopyran-2-one, a natural product found in many plants, has been scrutinized as a potentially potent cholinesterase inhibitor comparable to tacrine. Recent molecular docking simulations by Tanoli et al. reveal that the most potent coumarin derivatives, containing both piperidinyl and ethoxyl groups, lower the binding energy with AChE by nearly 1.5-fold that of tacrine. Large substituents attached to the coumarin ring enable these inhibitors to increase molecular contact with grooves in the enzyme active site gorge including, notably, simultaneous interactions with residues from the choline binding site, the peripheral aromatic site, and the catalytic triad [121]. Coumarin-linked thiourea derivatives exhibit similar binding modes in docking calculations, forming contact with the catalytic triads of both AChE and BChE [122]. However, in this same study, the thiourea group is observed to consistently hydrogen bond with Tyr146 of AChE, with no analogous hydrogen bonding observed for BChE due to the structural differences between the two enzymes. Moreover, hydrophobic interactions appear to dominate electrostatic interactions in these docking results. Another coumarin derivative, 7-hydroxycoumarin, also displays dual binding site capability with both the PAS and CAS of both cholinesterases [123]. A recent docking and MD study conducted by the Mubarek laboratory indicates that while hydrophobic interactions are dominant in stabilizing AChE in complex with coumarin derivatives, structural stability of BChE in complex with these species is predominantly due to hydrogen bonding [124].

#### 4.1.7. Other Pharmaceuticals

A variety of other FDA-approved drugs not intended to treat AD have been examined as possible cholinesterase inhibitors using computational methods. For example, adamantyl-based ester derivatives, which have been more widely used as acne, type 2 diabetes, and anti-viral medications, were studied in complex with both cholinesterases via docking; compounds with a methoxy substituent at position three on the phenyl ring showed the highest potential of binding strongly to both AChE and BChE [125].

Curiously, a number of marine metabolites, which show promising capabilities as pharmaceuticals, were docked with AChE and analyzed as potential treatments for AD [126]. From a database of FDA-approved drugs, Hassan et al. recently employed a screening technique that chose five drugs with higher capacity for AChE inhibition: Risperidone (for schizophrenia and bipolar disorder), Domperidone (for nausea and vomiting), Verapamil (for high blood pressure), Tamsulosin (for enlarged prostate), and Cinitapride (for nausea and ulcers); MD simulations indicated that all complexes were stable [127].

Another recent study by Ozer and coworkers employed docking calculations to understand interactions between BChE and fluoxetine, also known as Prozac, which is commonly used to treat anxiety, obsessive compulsive disorder, and anorexia [128]. Coupled with experimental kinetics measurements, fluoxetine proved to be a competitive inhibitor of BChE that binds deep in the active site gorge [128]. Previous MD and docking studies of BChE in complex with berberine derivatives, compounds found in medications for diabetes and high cholesterol, indicated that Trp82 (CBS), Gly117 (OAH), Trp231 (ABS), and Phe329 (CBS) are all important for binding [129]. Pyridoxine, commonly known as vitamin B6, is also recognized for its AChE inhibition ability: MD simulations of this complex indicate that the ligand forms a covalent bond with Ser203 of the catalytic triad, thus creating a steric barrier for acetylcholine [130].

#### *4.2. Narcotics*

Another field of inhibition centers around the interactions between narcotics and the cholinesterases, particularly BChE. As mentioned above, BChE is a promiscuous plasma enzyme, and can thus hydrolyze a variety of substrates during circulation. For example, BChE is one of the enzymes responsible for hydrolyzing, and subsequently activating, heroin, making it an ideal enzyme to target as a treatment for heroin overdose. In a recent study by Zhou et al. that addresses the need to block the activation of heroin by BChE, novel inhibitors from solanaceous alkaloid scaffolds were discovered via virtual screening, thus identifying a series of highly selective BChE inhibitors [131].

Recent docking and MD studies have also investigated cholinesterase interactions with nicotine and numerous derivatives thereof. For example, investigations into nicotine-AChE complexes using MD simulations, in tandem with experimental efforts, showed R-nicotine to more strongly disturb the secondary structure of, and to be a stronger inhibitor of, AChE than the S-nicotine analog [132]. Nicotine is the parent compound to the neonicotinoid family, which are present in commercial insecticides and, unsurprisingly, act as agonists to ACh receptors. In a 2018 study, Terali assessed the seven commercially available neonicotinoids using docking calculations, revealing different binding modes with AChE and BChE, and suggesting them to be potential compounds to treat cholinergic and non-cholinergic AD pathogenesis [133].

As discussed in the catalysis section above, cocaine is also a narcotic of interest for cholinesterase studies, and details of the catalytic mechanism of cocaine hydrolysis were discussed in that section. From MD simulations and hydrogen bonding energy (HBE) calculations, the energy barrier for hydrolysis of ACh by each enzyme was compared to the hydrolysis of (+)- and (−)-cocaine by human BChE and both of these differences were found to be approximately 3–5 kcal/mol [62]. This energy difference is attributed to the fact that only Gly117 and Ala199 in the OAH of BChE interacts with the carbonyl oxygen of cocaine, with Gly116 not participating [62,134].

A focal point in computational research of BChE-cocaine complexes is how BChE can be mutated to hydrolyze cocaine faster and remain in circulation longer as a possible treatment for cocaine overdose [30]. For example, MD simulations uncovered that mutations of Phe547, Met554, and Phe561 (in the C-terminus section of hBChE) to more hydrophobic residues may increase its circulation [30]. Circulation time of BChE may also be increased due to the introduction of more cross-subunit disulfide bonds, resulting in higher dimer stability as suggested in an MD study by Fang et al. [135].

Decreasing the activation energy for transition states, particularly for the rate-determining step, is another approach to amplify BChE catalytic efficiency. After combined computational and experimental studies revealed the rate determining step for the Ala328Trp/ Tyr332Ala and Ala328Trp/Tyr332Gly BChE mutants [78,136], MD simulations and virtual screening techniques were employed to determine how these mutations could lower these energy barriers, eventually yielding a mutant that was approximately 2000-fold more catalytically efficient than WT BChE, the 5-point mutant Ala199Ser/Phe227Ala/Ser287Gly/ Ala328Trp/Tyr332Gly [137]. The aforementioned Ala328Trp/Tyr332Ala and Ala328Trp/ Tyr332Gly BChE mutants were identified as being more catalytically active than WT human BChE, whose binding modes with (−)- and (+)-cocaine isomers as prereactive and nonprereactive complexes are displayed in Figure 5 [138]. The Ala328Trp/Tyr332Ala/Tyr419Ser mutant studied therein lost catalytic potency, as the conformation in which cocaine binds to this mutant is not conducive to catalysis [138]. Not all BChE mutants, however, are more catalytically efficient than WT BChE. Prompted by kinetics studies [139,140], a 2015 MD study of mouse and human BChE, alongside their Ala199Ser/Ser227Ala/Ser287Gly/Ala328Trp/ Tyr332Gly mutants, showed that the parent enzymes are approximately 250-fold more catalytically effective than their derivative hydrolases [141].

−

− − **Figure 5.** Wild-type BChE bound to (**a**) (−)-cocaine in a non-prereactive complex, (**b**) (−)-cocaine in a prereactive complex, (**c**) (+)-cocaine in a non-prereactive complex, and (**d**) (+)-cocaine prereactive complex [138].

#### *4.3. Organophosphates*

#### 4.3.1. Reversible Inhibition

Many organophosphates (OPs) with good leaving groups (weak bases) bind irreversibly to both cholinesterases, and molecular recognition complexes have been examined with docking and MD simulations. Human, mouse, and housefly models of AChE, as well as horse BChE, in complex with O,O-dialkylphosphate inhibitors suggest that the amino acid at position 400, which is either valine or phenylalanine, plays a key role in determining how well the OP will bind in the pocket based on the bulkiness of that amino acid [142]. Furthermore, Lee and Barron studied insect and mouse AChE with OP inhibitors using docking/QSAR techniques, demonstrating that Leu328 in the acyl binding site of insect AChE allows for less enzyme specificity in comparison to Phe295 of mouse AChE, and emphasized the importance of-interactions with the OAH [6].

Around this same time, Veselinovic et al. used Monte Carlo as part of their QSAR analysis to identify the best AChE inhibitors from a database of 278 OP compounds, with the goal of reducing cholinergic activity [143]. In recent months, Yang et al. published their combined Monte Carlo/MD study of *Tc*AChE adsorption in charged monolayers, which revealed that binding sites in the active site gorge orient themselves toward positively charged surfaces and away from negatively charged surfaces, a somewhat intuitive result given that cholinesterases have evolved to attract and hydrolyze positively charged choline moieties, but also providing useful insight for experimentalists using AChE as a means to detect OP compounds [144]. Also reported in the last few months were docking calculations of *Electrophorus electricus* AChE in complex with the voluminous and negatively charged 12-tungstosilicic acid and 12-tungstophosphoric acid, which allowed for detection of a previously-unknown allosteric binding site that has been subsequently labeled β-AS [145].

BChE-OP molecular recognition complexes have also been investigated by the Sorin laboratory using docking and MD methods. In a collaborative 2017 study that featured experimental work, the structural basis for relative *K*<sup>I</sup> values was probed via massive docking calculations for an assortment of dialkyl and aryl phosphate inhibitors in complex with BChE [146]. That same year, a massively-parallel MD study involving a very limited number of dialkyl phenyl phosphates probed the entropy change associated with binding to these OPs, and demonstrated there to be residual entropy associated with larger, more complex inhibitors that can sample from a much broader array of binding microstates (poses), thus adding to the stability of those larger and more flexible inhibitors entropically [147].

More recently, this same laboratory studied an array of dialkyl phenyl phosphate inhibitors via MD simulation, including numerous phenyl substitutions that had been previously probed via docking calculations [146] and a small set of alkyl-to-cholinyl substitutions to mimic the chemistry of the natural substrates. It was noted therein that larger R-groups increase van der Waals contact area between the enzyme and the ligand, with S-enantiomers apparently binding more strongly than their R analogs [35]. In an effort to characterize the observed modes of binding in the flexible BChE-OP complexes, contact tables such as that shown in Figure 6 were used to highlight specific interactions between portions of the inhibitor and specific binding sites and amino acids in the BChE active site gorge. Each row below the label rows at the top represents a binding mode, and every column is an amino acid residue known to participate in binding. Here, contacts are identified as chemical groups separated by 5 Å or less and the cell entries report which functional group dominates that interaction based on relative intermolecular interaction strength.

In a 2020 follow-up study, these BChE-OP were subject to massively-parallel MD simulation and then used as model systems around which to develop a methodology for accurately identifying binding modes from such rich data sets [148]. The resulting technique, brute force *k*-means clustering of surface-weighted interaction fingerprints (SWIFs), employs simple and intuitive statistical criteria to identify binding modes, and bypasses the heuristic nature of the *k*-means clustering algorithm. The contact table in Figure 6, taken from this most recent publication [148], demonstrates distinct binding β

modes for the OP inhibitor that binds most strongly to BChE of those so far studied by that laboratory and their collaborators, with *K*<sup>I</sup> = 1(±0.4) µM. μ

**Figure 6.** Contact table for the DIM5 inhibitor in BChE binding pocket. SWIFs were taken post 80 ns from one thousand 110 ns MD simulations [148].

#### 4.3.2. Irreversible Inhibition, Activation, and Reactivation

While the section above centered on reversible molecular recognition of OP inhibitors, "aged" cholinesterases are those that have undergone phosphorylation and have experienced significant structural change, thereby rendering them catalytically impotent. A survey of the energy landscape of the acyl pocket loop uncovered that the products of the reaction between AChE and diisopropyl fluorophosphate deviate significantly from the AChE crystal structure [149].

Like some narcotic and catalysis studies, researchers have focused on how BChE mutants can address and add insight to OP poisoning. For instance, Dwyer et al. found that BChE mutants Tyr332Ser, Asp340His, and Tyr332Ser/Asp340His all resist nerve agents by modifying the size of the "main door" [150]. However, in the case that the enzyme is already bound to the OP, Masson et al. focused on transition states in order to find more catalytically efficient mutants [151], and a more recent study by this same team showed that other mutations, such as Asn322Glu/Glu325Gly with an alternate Ser198, His438, Asn322Glu catalytic triad, allow complexes to self-reactivate [60].

MD simulations have also shown that minor mutations, such as replacing Gly116 in the OAH of BChE, could cause severe structural deviations [152], and simulations of the Gly117His and Gly117Asp mutants of BChE revealed the presence of a water molecule near Ser198 in the catalytic triad, which may be responsible for the reordering of water and subsequent conformational changes present in the mutants [152]. When these water molecules are replaced by glycerol molecules in cresyl saligenin phosphatephosphorylated BChE, there is a conformational change caused by His438 leading to a less reactive intermediate complex [5]. The structure of the Gly117His BChE mutant, studied via X-ray by Nachon et al. [153], was probed by Amitay and Shurki, who identified a single conformation from a set of computationally-generated structures that would fully reproduce the acetylation of acetylthiocholine [154]. The Gly117His BChE mutant was further explored as a potential OP bioscavenger in a QM/MM study by Yao et al. that reported improved activity against sarin by reducing the rate-determining energy barrier compared to wild type hBChE [155], as demonstrated in Figure 7.

The introduction of another compound as a reactivator to aged, or phosphorylated, cholinesterase has also been considered. Many reactivation studies have centered around AChE because of its crucial role in neurotransmission, which is interrupted by the introduction of nerve agents and insecticides, but a recent structure-based study suggested that small molecules (<200 Da), such as oximes, may act as strong antidotes for OPpoisoned AChE [9]. Oximes have been widely explored as reactivators for deactivated cholinesterases, so it is no surprise that the structure-activity relationship and docking studies proposed oxime-based compounds as antidotes for cholinesterases poisoned by OPs [156,157].

Nerve agents such as VX (venomous agent) and sarin (GB, as classified by the US-American military) have also been a focal point of studies on cholinesterase reactivation. Oximes used as reactivating compounds for OP-poisoned AChE are supported by computa-

tional and experimental research [158]. The prereactive complex of the oxime HI-6 and *Mus musculus* AChE covalently inhibited by sarin was examined by Allgardson et al., whose X-ray investigations and DFT calculations provided an essential foundation for research into the reactivation mechanism of OP-poisoned AChE [159]. Monte Carlo calculations by Veselinovic et al. of AChE-sarin reactivation reiterated that pyridinium oximes are decent antidotes [160]. A more recent study of charged and uncharged oximes by de Souza et al. compared these species with VX- and GB-poisoned AChE: while charged oximes proved to outperform the uncharged oximes, it is also an unfortunate reality that charged oximes do not cross the blood-brain barrier very well, making physical intake of the better reactivator more difficult [161]. Despite this setback, oximes are generally explored in more depth compared to pre-exposure antidote carbamates because carbamates also change the AChE structure via carbamylation [162].

**Figure 7.** Free energy profiles for the reactivation of sarin-phosphorylated human BChE for (**a**) wild type hBChE and (**b**) the G117H mutant discussed in the text [155].

Tabun (GA, as designated by the US-American military) is, unlike other nerve agents, particularly resistant to oxime compounds as reactivators [163]. This resistant quality has motivated researchers to find more effective oxime derivatives for tabun-cholinesterase complex reactivators. Dimethyl(pyridin-2-yl)sulfonium based oximes were examined at the DFT M05-2X/6-31G\* level and determined to be better reactivators, as they lower the energy barrier by 4.4 kcal/mol [164], and hierarchical ab initio calculations revealed that charged oxime derivatives as antidotes to tabun bound AChE are stronger than normal oxime compounds due to specific stereoelectronic characteristics [163]. Indeed, a 2014 study by Lo and Ganguly found charged oximes to be more effective than their uncharged analogs, and their QM/MM studies further suggested that N-(pyridin-2-yl)hydroxylamine is a better antidote than traditional oxime treatments and that it has a similar IC<sup>50</sup> value [165].

Treatments for general nerve agent and insecticide poisoning have utilized oxime derivatives as well. Reactivation of a VX-AChE complex using a deprotonated pralidoxime, or 2-pralidoxime (2-PAM), occurs through consecutive addition-elimination steps and shows promising results as an antidote [166]. Docking and QM/MM methods paired with experimental observations revealed that trimedoximes show potential to reactivate *Mus musculus* AChE, with the AChE-VX complex showing the best results [167], and MD simulations of 2-PAM with phosphorylated AChE support this claim [168]. The importance of protonated Glu202 in the reactivation of VX-inhibited mouse AChE was observed in QM/MM simulations performed by Driant et al. [169]. Further, symmetrical and unsymmetrical isoquinolinium-5-carbaldoximes showed strong inhibition for both cholinesterases; the weaker inhibitors were selected for additional experimental and computational investigation [170]. Interestingly, QSAR studies found that a combination compound consisting of tacrine and aroylacrylic acid phenylamide moieties showed potential as pre-exposure OP-poisoning antidotes [171].

#### *4.4. Other Organic Moieties*

#### 4.4.1. Hydrocarbons

The Sepˇci´c laboratory studied the interactions of the carbon-based nanomaterials (NM) carbon black (CB), fullerene (C60), and graphene oxide (GO) in complex with AChE experimentally and with docking and MD simulations, finding that CB inhibited AChE most efficiently, while C<sup>60</sup> was least efficient and interactions with the GO surface allowed AChE to retain its native shape and activity [172].

Flavonoids are targeted as potential inhibitors that are not regulated by the FDA. Vats et al. found a number of flavonoid analogues to be novel AChE inhibitors via QSAR analysis [173]. Another sub-class of flavonoids are catechins, including hydroxyl-rich epicatechin, which has been undergoing trials as a potential therapeutic for diabetes and cancer. Of these, epicatechin 3,5-O-digallate was investigated with docking and MD simulations in complex with BChE and found to bind closely to the His484 residue of the catalytic triad with as many as six stabilizing hydrogen bonds [174].

As noted above, investigations into cholinesterase structure and function in the presence of certain toxins are of significant interest. The behavior of aflatoxin, for example, which is regulated by the FDA, not as a pharmaceutical but as a toxin, was examined by Sanson et al. in complex with AChE using MD simulations [50], which revealed that the presence of aflatoxin and its interaction with Trp84 (in the CBS) caused enlargement of the active site gorge [50]. Furthermore, a recent study by the de Almeida laboratory showed aflatoxin M1, a toxic natural compound found in contaminated dairy products, to be another potential inhibitor of AChE, which binds to the CAS region of AChE but does not bind to or inhibitor BChE [175].

Other compounds that have been considered as BChE inhibitors include derivatives of 2-phenylbenzofuran [176]. MD simulations revealed that, while these derivatives bound to the PAS and CAS sites quite well, a derivative that included a para-position hydroxy group on the phenyl moiety improved inhibition against BChE [176]. Phenyl valerate is another aromatic hydrocarbon that has been studied recently in complex with BChE by Estevez et al.: they observed via MD simulation that phenyl valerate inhibits BChE at different ends of the active site gorge, thereby inhibiting the hydrolysis of ACh; it was experimentally determined, however, that both phenyl valerate and ACh can be hydrolyzed simultaneously [177], highlighting the need for extensive simulation time with respect to complex systems and caution when interpreting the results of those simulations.

#### 4.4.2. Nitrogenous Compounds

#### Amines

Amines are one of the most common functional groups occurring in nature. As a reminder, compounds are placed in this section due to the general abundance or consistent occurrence of amine groups in the inhibitor series studied in a given paper. For example, piperidine is an amine-containing compound that has been considered as a cholinesterase inhibitor in a handful of publications, including phenoxyethyl piperidine derivatives, which were studied via MD simulation in complex with electric eel AChE and horse BChE [178]. The phenoxyethyl derivatives most structurally similar to donepezil had the ability to bind to both the CAS and PAS, while many others bound only to the CAS [178]. In contrast, piperidine compounds substituted with arylaminopropanone were recently examined by Hudcova et al. via docking, MD, and QM/MM approaches and were compared to rivastigmine and galantamine [179]. In fact, our current appreciation for the role of hydrophobic active sites residues in binding is highlighted by a previous study by Khayamian and coworkers, who tested 68 piperidine and amine compounds as AChE inhibitors using docking and MD simulations, thereby confirming that hydrophobic interactions are a dominant factor in cholinesterase-inhibitor binding [180].

A handful of other amine compounds have been considered as possible cholinesterase inhibitors. Docking simulations of 4-acetoxy-plakinamine B in complex with AChE uncovered that the inhibitor binds primarily with the PAS and ABS [181]. Moreover, Shrivastava

et al. recently studied 23 p-aminobenzoic acid derivatives in complex with both AChE and BChE, which were compared to the binding affinity of donepezil [182], and 1Hbenzimidazole compounds with amine substituents were shown, in a previous study, to prefer BChE in docking and MD simulations [183]. Indeed, around 85 amine-containing compounds were used as an input to a QSAR study by Abuhamdah et al., with 24 compounds exhibiting micromolar IC<sup>50</sup> values [184].

#### Amides, Imides, Imines, and Carbamates

Other miscellaneous but common nitrogenous functional groups and compound types are amides, imines, and carbamates. Imides are one of the least examined groups, which is acceptable given the large swath of papers with inhibitor series that have little in common with each other. A QSAR analysis of 84 N-aryl-monosubstituted derivates provided 42 imide inhibitors and emphasized the importance of-interactions with the Trp82 and Trp86 residues in BChE and AChE, respectively, via MD simulations [185]. Imines are close behind, with N-(1-(5-bromo-2-hydroxyphenyl)-ethylidene)-3,4,5-trihydroxybenzohydrazide, a Schiff base derivative and the most potent AChE inhibitor in this series of compounds, revealed via docking calculations to interact mostly with the PAS and ABS [186], as was reported above.

Amide compounds have been more widely studied as cholinesterase inhibitors, including anandamides and acylethanolamides (NAEs), commonly found in most tissues, along with oleoylethanolamide and palmitoylethanolamide, which were docked with BChE; the latter were found to be uncompetitive inhibitors of BChE, and anandamides were found to be noncompetitive [187]. Docking and MD were also utilized to study the inhibitory activity of 4-aryl-oxo-2-aminylbutamides with both AChE and BChE. Although many of these compounds were ineffective towards AChE, the most potent AChE inhibitors displayed a tendency for stronger interactions between -NH moiety and the Tyr124 hydroxyl in the PAS [188]. More recently, Singh and Gupta used a multi-pronged QSAR analysis, docking, and MD approach to study a series of potential AChE inhibitors in which over half contained amide groups, reporting that inhibitors lacking amide groups received lower docking scores [189].

Carbamates are structurally similar to amides and carbamate-based inhibitors thus behave similarly to amide-based inhibitors. Recently reported RMSD calculations of carbamate-based inhibitors in complex with AChE were found to reach structural equilibrium after about 6 ns of MD simulation time and demonstrated that carbamate inhibitors with aromatic rings were more strongly drawn to AChE's binding pocket [190]. Indeed, from another recent MD-based study of thymol carbamates in complex with BChE, the importance of hydrophobic interactions with Trp82 (CBS), Gly116 (OAH), and Gly197 was emphasized in relation to MR-complex stability, as were water mediated interactions within the complex [191]. Analogously, the stability of AChE-inhibitor complexes was found to depend strongly on hydrophobic interactions with Tyr341 and Trp286 (near gorge entrance) and hydrogen bonds with Tyr124 (also near gorge entrance) and Phe295 (farther down inside the gorge) [191].

#### Nitrogenous Heterocyclic Rings and Derivatives

One class of molecule common to this group are indoles and indole derivatives, a number of which have been studied via docking. For example, docking calculations by Dominguez et al. with AChE revealed that meta-substituted benzylamine indole derivatives outperformed other indole cholinesterase inhibitors [192]. Dileep et al. also docked indole-3-acetic acid (IAA) and indole 3-butyric acid (IBA) derivatives with AChE; these indoles are known as auxins, or plant growth regulators, which are used in culture experiments for plant tissue [193]. Most recently, Bingul et al. docked six 4,6-dimethoxyindole based hydrazide-hydrazones that were found to bind to both AChE and BChE more strongly if they contained a phenyl group [194].

Piperine, a chemical closely associated with black pepper, has also been investigated as a potential cholinesterase inhibitor. Arylaminopropanone derivatives substituted with piperidine were docked and simulated with both cholinesterases, proving that they perform most similarly to galantamine and rivastigmine among the series of arylaminopropanone with N-phenylcarbamate moieties [179]. In another recent docking study, piperidine and curcumin, a chemical found in turmeric plants, were found to bind most closely to AChE, with binding energies of −10.5 kcal/mol and −9.6 kcal/mol, respectively [195], and it has been shown that cholinesterase binding sites most responsible for strong interactions with piperidine, as with many of the species discussed above, are the CAS and the PAS regions of the gorge [196,197].

Docking and MD simulations revealed that isoalloxazine derivatives with an ortho dimethoxybenzyl group were favorably bound to the peripheral binding site of AChE [198]. AChE-docked 4-aminopyridine semicarbazone derivatives with biphenyl rings had stronger hydrophobic interactions and overall better binding affinities than derivatives without those groups [199]. Carbazole-based stilbene derivatives were also docked with both cholinesterases, as well as the Aβ1-42 peptide, and displayed potential as a multitarget inhibitor for AD in a 2020 study by Patel et al. [200].

One interesting potential cholinesterase inhibitor comes from the Chilean *Rhodophilia* (Amaryllidaceae) plant. Docking simulations were performed on *Rhodophilia* compounds with the highest alkaloid compositions, with IC<sup>50</sup> values also reported from in vitro experiments [201]. The Bohorquez laboratory also studied novel N-allyl/propargyl 4-substituted 1,2,3,4-tetrahydroquinoline derivatives in complex with both AChE and BChE. Their results from docking, MM/GBSA simulations, and experimental work showed a high correlation between the calculated binding free energy and inhibitor activity for both cholinesterase targets [202]. Several years later, this same team studied tetrahydroquinoline(THQ) isoxazole/isoxazoline compounds, which proved to have similar binding modes to galantamine [203]. Furthermore, tricyclic and tetracyclic nitrogen bridgehead compounds in complex with AChE were investigated with docking and MD simulations as potential inhibitors for AD treatment by the Decker group, with their strongest inhibitor in the tens of nanomolar regime [204].

Interestingly, numerous molecules are studied in tandem with other receptors and enzymes. For example, the Decker group effort noted above also used docking and MD to understand the high affinity of their inhibitors for the human histamine H<sup>3</sup> receptors [204] and Samadi et al. docked heterocyclic substituted alkyl and cycloalkyl propargyl amines to AChE, BChE, and monoamine oxidases, which are also closely linked to AD [205]. Similarly, small benzimidazole-based molecules were studied as both BChE inhibitors and as human cannabinoid receptor agonists, highlighting the possible future development of a dual-acting therapeutic to treat AD [206].

#### 4.4.3. Organosulfates

It is noteworthy that many studies have used inhibitor series that consistently contain sulfur groups. In some cases, sulfur groups are not the intended focal point; however, since sulfur is a common element in nature, it is deserving of a sub-section. One example is a pair of recent studies by Hassan et al., who studied a series of amide and piperazine sulfonamide derivatives via docking with AChE and BChE [207,208], in which the piperazine sulfonamide derivatives with substituted alkenes proved to be the strongest BChE inhibitors [208]. Another recent example comes from the 2019 study of Yang et al., who looked at a series of inhibitors discovered from structure-based pharmacophore virtual screening and included a number of potential inhibitors containing either thiols or sulfones [209].

Phthalimide-dithiocarbamate hybrids were also recently docked with both cholinesterases, revealing binding modes that are comparable to the FDA-regulated donepezil and rivastigmine pharmaceuticals [210]. In addition, an earlier study showed 7H-thiazolo[3,2-b]-1,2,4 triazin-7-one derivatives with two substituents on the phenyl group closest to the sulfur to show promise of dual CAS and PAS binding and AChE inhibition [211]. Interestingly, a pattern of correlation was found between decreasing fluorescence intensity and increasing binding activity of thioflavin-T with AChE, which was confirmed with docking and MD simulations [212].

Clearly, compounds with thiol or other sulfur-containing groups can be quite effective cholinesterase inhibitors, and early MD simulations revealed that a rivastigmine analog with a sulfur system was 192-fold more efficient at inhibiting cholinesterases reversibly than the rivastigmine parent molecule [113]. MD simulations also demonstrated that the non-competitive substrate acetylthiocholine inhibits both cholinesterases at different ends of the active site; despite the partial competition, experimental work points out that acetylthiocholine and ACh can be hydrolyzed virtually simultaneously [177]. The last compound we will mention here, benzothiazepine, was preferentially bound to BChE due to stronger hydrogen bonding, as observed in MD simulations [11].

#### *4.5. Proteins, Nucleic Acids, and Salts*

#### 4.5.1. Protein and RNA Binding

Cholinesterase inhibition by small proteins and RNAs has also been studied computationally for a number of species. For instance, fasciculin II is a peptidic three-finger snake toxin and was observed in 5 ns MD simulations to bind to the mouth of the active site gorge in mouse AChE, effectively blocking any substrate from entering the gorge [213], as was observed in an X-ray structure from the Sussman laboratory just a few years earlier [214]. AChE was also modeled with cytochrome c (Cyt c), a hemeprotein generally associated with respiratory cell functions, and it was reported that AChE interactions with Cyt c play a crucial role in apoptosome formation. Macro-modeling studies reveal that Cyt c binds to the PAS of AChE, blocking gorge access a la fasciculin, and also that binding modes with AChE are similar regardless of whether the heme group in Cyt c is present (Holo) or absent (Apo) [215].

Sohail and Rashid used docking and MD simulation to study interactions between the RNA recognition motif (RRM), the most abundant RNA-binding protein domain, and BChE. While it is unclear the degree to which RRM-bound BChE would be catalytically efficient, the authors note that gaining a better understanding of these interactions in vivo could prove highly useful in therapeutic development [216]. A recent docking-based follow-up to that article reports the use of microRNA (miR-132) as a potential inhibitor of AChE, with miR-132 binding predominantly via interactions with the catalytic triad, and thus blocking access to the substrate [217].

#### 4.5.2. Nucleobase Derivatives

Most nucleobase derivatives have focused on pyrimidine, a six-membered heterocyclic compound found in DNA and RNA. Examples of these derivatives include a series of di-phenylpyrimidine derivatives that were recently investigated as AChE inhibitors with docking and MD simulations [218,219]. The molecule labeled VB8 by Kumar et al. showed the highest activity against AChE by demonstrating additional (substituent-induced) interactions with the active site gorge [219]. Another recent docking and MD study reveals that uracil derivatives can inhibit both cholinesterases [220], and 6-methyluracil has also been recently modeled in complex with AChE and BChE, showing stronger binding than donepezil and stabilizing secondary binding to the PAS [7,221].

#### 4.5.3. Ion and Salt Binding

Although not all salts are organic, it would be irresponsible to exclude them from this review. A 2020 docking study by Yigit et al. found amine-tethered benzimidazolium salts with a trimethyl benzyl ring to be efficient AChE inhibitors due to the close proximity and strong interaction with the CAS and the PAS [222]. Previously, 2-*N*,*N*dimethylaminecyclohexyl 1-*N*′ ,*N*′ -dimethylcarbamate isomers and their methylsulfate salts were investigated computationally, in tandem with experimental work, as cholinesterase inhibitors, revealing that the lowest binding rate was 55% and the highest binding rate was 90% with BChE [223].

#### **5. Virtual Screening**

After exploring the multitudes of cholinesterase inhibitors that have been studied via computation, it is important to recognize progress and state-of-the-art improvements that have been made in the area of virtual screening, particularly as developed for or demonstrated on AChE and BChE. For those not familiar with the term, virtual screening is a computational method by which large libraries of small molecules can be searched for potential matches to a specific drug target. Many virtual screening methods that have been developed specifically in relation to the cholinesterase enzymes. As has been highlighted in numerous sections above, virtual screening studies are typically QSAR studies paired with docking and/or MD simulation; such studies are designed to find chemical entities and assess their potential inhibitory effectiveness [224].

As an example, Discovery Studio 2.5.5. was used to construct pharmacophore models with the goal of finding molecules that inhibit AChE and protect it from amyloid beta toxicity: from a sample of 62 compounds, only nine were found to interact favorably with AChE [225]. Some years later, the same research team used virtual screening to identify BChE-specific inhibitors from commercial databases of comprising 3.9 million compounds: virtual screening, docking, and bioassay reduced the list of possible matches to just six compounds [226]. Similar results have been reported for other virtual screening studies, including a computer-aided workflow that utilized hierarchical, structure-based screening, thereby yielding five potential cholinesterase inhibitors [227]. Furthermore, six inhibitors, three each for AChE and BChE, out of four commercial compound databases were found using structure-based pharmacophore modeling intended for AD treatment [209]. Similarly, both ZINC (zinc15.docking.org) and DrugBank (go.drugbank.com) were screened for reactivation oxime compounds [228]. As a true success story, the Gobec group developed a successful virtual screening method for BChE [229] that recently helped to realize the discovery of some of the most powerful reversible inhibitors of BChE known, with inhibitions constants in the picomolar to nanomolar range [230].

There are indications, however, that advancements in virtual screening remain to be made. For example, an early report of machine learning as a guide to virtual screening produced cholinesterase inhibitors with a wide range of IC<sup>50</sup> values, suggesting that certain proposed inhibitors may be too toxic [231]. Detecting false positives would also improve virtual screening studies, an issue that was addressed when screening compounds from Maybridge.com and ChemBridge.com databases for potential BChE inhibitors, in which five ligands were chosen after ADMET calculations, Lipinski's Rule of Five, and docking simulations [232]. Accounting for inhibitor stereochemistry with respect to active site geometry is another factor to be addressed. For example, virtual screening studies of 24 chiral organophosphates revealed that S-isomers exhibited stronger inhibitory activity towards AChE than their respective R-isomers [233].

Given their pronounced gorges and well-defined binding sites, along with the extensive literature on cholinesterase studies, it is no surprise that AChE and BChE have been used as example enzymes in a number of virtual screening studies. Some methods take hands-on approaches, such as utilization of a Monte Carlo approach paired with CORAL calculations [234]. Employment of steered MD simulations that calculate the work needed to remove the ligand from the binding site proved to be more efficient than conventional MD simulations and subsequent calculations of binding energies [235]. Using docking simulations as a means of assessing ligand mobility, a factor that is considered in virtual screening, was also demonstrated on BChE [236]. Automated docking software, such as ICM-Pro [146] and AutoDock [237], has also proven insightful in recent studies.

In contrast, other researchers have developed their own virtual screening algorithms to focus solely on the cholinesterases. For example, SHAFTS (SHApe FeaTure Similarity) is a 3D similarity calculation designed for AChE ligand discovery [238] and LiSiCa (Ligand Similarity using Clique algorithm) is a virtual screening development featuring BChE [239]. Lastly, the ADAM&EVE virtual screening method was presented with a focus on AChE inhibition and identified thirteen potential compounds from an original database

of 160,000 [240]. It will indeed be exciting to see the directions in which this area grows in the years to come.

#### **6. Conclusions**

Cholinesterase structure, function, and inhibition have proven to be a source of great interest for the application of a broad spectrum of modeling and computational approaches. This review examined articles that detailed the mechanisms by which substrates and inhibitors locate and enter the gorges of AChE and BChE, how specific binding sites within the active site gorge of these enzymes respond to and interact with specific gorge binding sites, and pathways that the products of hydrolysis and other small molecules may find to enter or leave the active site. Also reviewed were computational studies of the mechanisms and thermodynamics of substrate hydrolysis by both enzymes, revealing that AChE and BChE have distinct rate-determining steps. Catalysis studies directed at BChE hydrolysis of cocaine, and the various mutations that can speed up the catalytic process, were also explored in depth.

The bulk of our findings, however, revolved around inhibition. Some major pharmaceutical compounds such as tacrine, galantamine, donepezil, rivastigmine, quinazoline, and coumarin—most of which are conventional, well-known medications for the treatment of AD and other human ailments—and their derivatives, were docked or simulated with AChE and BChE. A sizeable group of additional FDA-regulated compounds were also discussed, with the narcotic section highlighting findings regarding the cholinesterases either bound to or hydrolyzing heroin, nicotine, and cocaine. Organophosphates were heavily explored above, including: discussion of various OP compounds as potential reversible inhibitors; the irreversible binding of nerve agents such as sarin, tabun, and VX; and enzyme reactivation by molecules such as oxime derivatives. The final inhibition section, and by far the most diverse, included other organic compounds including various hydrocarbons, and nitrogenous compounds such as amines, amides, carbamates, and nitrogen heterocycles. Other organic molecules discussed were organosulfates, protein and nucleic acid derivatives, and ionic inhibitors.

Finally, advancements in virtual screening methodologies and software were discussed, with some of these methodologies specific to cholinesterases and others simply featuring them as prime applications. Although there have been hundreds of computational cholinesterase studies published in the last two decades, it is exciting to consider the many possible directions that computation will lend itself to improve our understanding of these enzymes and their function in the future, with potential real-world applications in human disease therapies, treatments for pesticide and nerve agent poisoning, and management of drug overdoses.

**Author Contributions:** Conceptualization, D.D.B., N.N., and E.J.S.; investigation, D.D.B., N.N., J.M. (Jia Mao), and J.M. (Jessica Moore); writing—original draft preparation, D.D.B. and N.N.; writing review and editing, D.D.B., N.N., and E.J.S.; visualization, N.N. and E.J.S.; supervision, E.J.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This Review was supported by funding from the National Institute of General Medical Sciences of the National Institutes of Health under Award Numbers UL1GM118979, TL4GM118980, and RL5GM118978. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** N.N. and J.M. (Jia Mao) are grateful for undergraduate research scholarships from Women & Philanthropy. N.N. acknowledges a Boeing scholarship, and N.N. and J.M. (Jessica Moore) thank the CSULB College of Engineering for scholarship support. D.D.B., J.M. (Jessica Moore), and E.J.S. acknowledge support from the NIH BUILD program: this review was supported by funding from the National Institute of General Medical Sciences of the National Institutes of Health under Award Numbers UL1GM118979, TL4GM118980, and RL5GM118978. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Molecular Modeling Studies on the Multistep Reactivation Process of Organophosphate-Inhibited Acetylcholinesterase and Butyrylcholinesterase**

**Jakub Jo ´nczyk <sup>1</sup> , J ˛edrzej Kukułowicz <sup>1</sup> , Kamil Ł ˛atka <sup>1</sup> , Barbara Malawska <sup>1</sup> , Young-Sik Jung 2,3 , Kamil Musilek 4,5,\* and Marek Bajda 1,\***


**Abstract:** Poisoning with organophosphorus compounds used as pesticides or misused as chemical weapons remains a serious threat to human health and life. Their toxic effects result from irreversible blockade of the enzymes acetylcholinesterase and butyrylcholinesterase, which causes overstimulation of the cholinergic system and often leads to serious injury or death. Treatment of organophosphorus poisoning involves, among other strategies, the administration of oxime compounds. Oximes reactivate cholinesterases by breaking the covalent bond between the serine residue from the enzyme active site and the phosphorus atom of the organophosphorus compound. Although the general mechanism of reactivation has been known for years, the exact molecular aspects determining the efficiency and selectivity of individual oximes are still not clear. This hinders the development of new active compounds. In our research, using relatively simple and widely available molecular docking methods, we investigated the reactivation of acetyl- and butyrylcholinesterase blocked by sarin and tabun. For the selected oximes, their binding modes at each step of the reactivation process were identified. Amino acids essential for effective reactivation and those responsible for the selectivity of individual oximes against inhibited acetyl- and butyrylcholinesterase were identified. This research broadens the knowledge about cholinesterase reactivation and demonstrates the usefulness of molecular docking in the study of this process. The presented observations and methods can be used in the future to support the search for new effective reactivators.

**Keywords:** molecular modeling; reactivators; reactivation process; organophosphates; docking studies; acetylcholinesterase; butyrylcholinesterase

#### **1. Introduction**

Organophosphorus compounds (OPs) have been known for years as effective pesticides used in agriculture (e.g., paraoxon and chlorpyriphos) and as chemical weapons (e.g., sarin, venomous agent X (VX), tabun, and soman) belonging to the so-called "nerve agent" class representing some of the strongest man-made poisons [1–3]. The toxic effects of OPs result from irreversible blockade of acetylcholinesterase (AChE), the enzyme responsible for the degradation of acetylcholine [1,3]. The accumulation of this neurotransmitter in

**Citation:** Jo ´nczyk, J.; Kukułowicz, J.; Ł ˛atka, K.; Malawska, B.; Jung, Y.-S.; Musilek, K.; Bajda, M. Molecular Modeling Studies on the Multistep Reactivation Process of Organophosphate-Inhibited Acetylcholinesterase and Butyrylcholinesterase. *Biomolecules* **2021**, *11*, 169. https://doi.org/ 10.3390/biom11020169

Academic Editor: Jan Korabecny Received: 21 December 2020 Accepted: 22 January 2021 Published: 27 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

reactivators [17,18].

*Biomolecules* **2021**, *11*, x FOR PEER REVIEW 2 of 21

of OPs result from irreversible blockade of acetylcholinesterase (AChE), the enzyme responsible for the degradation of acetylcholine [1,3]. The accumulation of this neurotransmitter in cholinergic synapses within the central and peripheral nervous system causes, among other effects, muscle tremors, spasms, cardiac dysfunction, and respiratory failure, often leading to death [4,5]. In addition to blocking AChE, OPs also have the ability to inhibit butyrylcholinesterase (BuChE), which also degrades other choline esters in addition to acetylcholine. The World Health Organization (WHO) estimates that approximately 3 million people are accidentally or intentionally intoxicated each year by OPs used as pesticides, with approximately 220,000 of them being poisoned fatally [6]. Although usage of chemical weapons is internationally banned, there are unfortunately still cases of the deliberate use of OPs, such as the terrorist attacks in the Tokyo subway in 1995, the bombing of civilians in Syria from 2013 to date, the use of a Novichok agent

9]. Treatment options for OP intoxication are limited. They are based on the administra-

cholinergic synapses within the central and peripheral nervous system causes, among other effects, muscle tremors, spasms, cardiac dysfunction, and respiratory failure, often leading to death [4,5]. In addition to blocking AChE, OPs also have the ability to inhibit butyrylcholinesterase (BuChE), which also degrades other choline esters in addition to acetylcholine. The World Health Organization (WHO) estimates that approximately 3 million people are accidentally or intentionally intoxicated each year by OPs used as pesticides, with approximately 220,000 of them being poisoned fatally [6]. Although usage of chemical weapons is internationally banned, there are unfortunately still cases of the deliberate use of OPs, such as the terrorist attacks in the Tokyo subway in 1995, the bombing of civilians in Syria from 2013 to date, the use of a Novichok agent against Sergei and Yulia Skripal in 2018, and the recent poisoning of Alexei Navalny [7–9]. Treatment options for OP intoxication are limited. They are based on the administration of atropine, which reduces the overstimulation of muscarinic receptors, benzodiazepines that reduce convulsions, and oximes, which are specific antidotes capable of reactivating AChE [5,10]. The effectiveness of oximes depends on the type of OP and the time elapsed from poisoning. The small number of such compounds currently available for use as drugs and their often low activity result in the need to search for new reactivators. Intravenous administration of BuChE, alone or in combination with its reactivator, as a pseudocatalytic bioscavenger of OPs in the bloodstream is being considered for use as a supportive treatment that allows more rapid clearance of OPs from the patient [11,12]. However, the lack of effective BuChE reactivators currently limits the implementation of this therapeutic solution. tion of atropine, which reduces the overstimulation of muscarinic receptors, benzodiazepines that reduce convulsions, and oximes, which are specific antidotes capable of reactivating AChE [5,10]. The effectiveness of oximes depends on the type of OP and the time elapsed from poisoning. The small number of such compounds currently available for use as drugs and their often low activity result in the need to search for new reactivators. Intravenous administration of BuChE, alone or in combination with its reactivator, as a pseudocatalytic bioscavenger of OPs in the bloodstream is being considered for use as a supportive treatment that allows more rapid clearance of OPs from the patient [11,12]. However, the lack of effective BuChE reactivators currently limits the implementation of this therapeutic solution. At the molecular level, the irreversible blockade of cholinesterases by OPs is based on the phosphorylation of serine, which is a part of the catalytic triad, also called the esteratic site, of the enzymes. This prevents the binding of the ester group of acetylcholine at this site and its further hydrolysis [1,13]. The other residues included in the catalytic triad are histidine and glutamic acid. The OP-enzyme complex may undergo an aging process. This process involves the dealkylation of an alkoxy moiety from the phosphorus

At the molecular level, the irreversible blockade of cholinesterases by OPs is based on the phosphorylation of serine, which is a part of the catalytic triad, also called the esteratic site, of the enzymes. This prevents the binding of the ester group of acetylcholine at this site and its further hydrolysis [1,13]. The other residues included in the catalytic triad are histidine and glutamic acid. The OP-enzyme complex may undergo an aging process. This process involves the dealkylation of an alkoxy moiety from the phosphorus atom of the OP. This process is clinically unfavorable because the aged complex is not susceptible to reactivation [14–16]. Other important structural elements of cholinesterases are the anionic site, which is responsible for binding the substrate quaternary ammonium group, and the acyl pocket, which binds the alkyl chains of acyl moieties. There is also an oxyanion hole that stabilizes the transition complex during the enzymatic reaction and a peripheral anionic site (PAS) located at the entry to the active gorge [17]. Notably, the active site in BuChE is significantly larger than that in AChE. This is due to the replacement of some of the aromatic amino acids within the acyl pocket, anionic site, and PAS in AChE with smaller aliphatic or even polar residues in BuChE. This translates into different substrate specificities of these enzymes as well as selective binding of inhibitors or reactivators [17,18]. atom of the OP. This process is clinically unfavorable because the aged complex is not susceptible to reactivation [14–16]. Other important structural elements of cholinesterases are the anionic site, which is responsible for binding the substrate quaternary ammonium group, and the acyl pocket, which binds the alkyl chains of acyl moieties. There is also an oxyanion hole that stabilizes the transition complex during the enzymatic reaction and a peripheral anionic site (PAS) located at the entry to the active gorge [17]. Notably, the active site in BuChE is significantly larger than that in AChE. This is due to the replacement of some of the aromatic amino acids within the acyl pocket, anionic site, and PAS in AChE with smaller aliphatic or even polar residues in BuChE. This translates into different substrate specificities of these enzymes as well as selective binding of inhibitors or The exact mechanism of cholinesterase reactivation by an oxime involves a nucleophilic attack of the oxime ion on the OP phosphorus atom bound to serine (SER) in the

The exact mechanism of cholinesterase reactivation by an oxime involves a nucleophilic attack of the oxime ion on the OP phosphorus atom bound to serine (SER) in the enzyme active site (Figure 1). As a result, a transient oxime-OP-SER adduct with a trigonal bipyramidal geometry is created. Subsequently, the OP-SER bond is broken, which restores the catalytic activity of cholinesterase and releases the OP-oxime complex. This complex should have low affinity for the binding enzyme pocket to quickly dissociate from it without rephosphorylation or further enzyme blockade [16,19,20]. enzyme active site (Figure 1). As a result, a transient oxime-OP-SER adduct with a trigonal bipyramidal geometry is created. Subsequently, the OP-SER bond is broken, which restores the catalytic activity of cholinesterase and releases the OP-oxime complex. This complex should have low affinity for the binding enzyme pocket to quickly dissociate from it without rephosphorylation or further enzyme blockade [16,19,20].

**Figure 1.** Mechanism of cholinesterase reactivation by oxime compounds. **Figure 1.** Mechanism of cholinesterase reactivation by oxime compounds.

Effective reactivators, in addition to a strong nucleophilic moiety such as the oxime moiety, must also have elements that allow interaction with the enzyme active site and

facilitate favorable positioning of the oxime fragment with respect to the OP phosphorus atom. These criteria are fulfilled by the strongest reactivators currently known, which are quaternary pyridinium cations with an aldoxime moiety substituted in position 2 or 4 of the pyridine ring. Due to the number of aromatic rings in the structure, they can be divided into monopyridinium and bispyridinium derivatives [21]. Examples of such reactivators are presented in Figure 2 [21–23]. facilitate favorable positioning of the oxime fragment with respect to the OP phosphorus atom. These criteria are fulfilled by the strongest reactivators currently known, which are quaternary pyridinium cations with an aldoxime moiety substituted in position 2 or 4 of the pyridine ring. Due to the number of aromatic rings in the structure, they can be divided into monopyridinium and bispyridinium derivatives [21]. Examples of such reactivators are presented in Figure 2 [21–23].

Effective reactivators, in addition to a strong nucleophilic moiety such as the oxime moiety, must also have elements that allow interaction with the enzyme active site and

*Biomolecules* **2021**, *11*, x FOR PEER REVIEW 3 of 21

**Figure 2.** Structures of selected examples of cholinesterase reactivators. **Figure 2.** Structures of selected examples of cholinesterase reactivators.

Although the crystal structures of both acetyl- and butyrylcholinesterase have been known for years, only some of them present a reactivator bound to the OP-enzyme complex. In addition, in almost all cases, the ligands have conformations in which the oxime moieties are located in the opposite direction from the OP fragment or are arranged in a geometrically unfavorable position that prevents nucleophilic attack on the OP phosphorus atom [24]. Allgardsson et al. [20] named this conformation nonreactive. They also obtained one crystal structure (PDB code: 5FPP) with the prereactivation conformation of the HI6 reactivator in chain A of the dimeric AChE complex inhibited by sarin. This structure presents the oxime moiety placed close to the phosphorus atom in a geometrically beneficial axial position. Interestingly, chain B obtained under the same conditions shows the nonreactivation conformation of asoxime (HI-6). It seems that the previous failures in the design of new effective reactivators may be partly caused by focusing only on the crystal structures that present reactivators in nonreactive conformations [25,26]. It is necessary to consider that an effective reactivator must be able to approach the phosphorus atom with the appropriate geometric constraints as presented by Allgardsson et al. [20]. Capturing this state by molecular modeling is not easy. It often requires the use of molecular dynamics techniques in addition to quantum methods because only they can reproduce the moment of oxime-OP bond formation [27–30]. Both molecular dynamics methods and, to an even greater extent, quantum calculations require high computational ca-Although the crystal structures of both acetyl- and butyrylcholinesterase have been known for years, only some of them present a reactivator bound to the OP-enzyme complex. In addition, in almost all cases, the ligands have conformations in which the oxime moieties are located in the opposite direction from the OP fragment or are arranged in a geometrically unfavorable position that prevents nucleophilic attack on the OP phosphorus atom [24]. Allgardsson et al. [20] named this conformation nonreactive. They also obtained one crystal structure (PDB code: 5FPP) with the prereactivation conformation of the HI6 reactivator in chain A of the dimeric AChE complex inhibited by sarin. This structure presents the oxime moiety placed close to the phosphorus atom in a geometrically beneficial axial position. Interestingly, chain B obtained under the same conditions shows the nonreactivation conformation of asoxime (HI-6). It seems that the previous failures in the design of new effective reactivators may be partly caused by focusing only on the crystal structures that present reactivators in nonreactive conformations [25,26]. It is necessary to consider that an effective reactivator must be able to approach the phosphorus atom with the appropriate geometric constraints as presented by Allgardsson et al. [20]. Capturing this state by molecular modeling is not easy. It often requires the use of molecular dynamics techniques in addition to quantum methods because only they can reproduce the moment of oxime-OP bond formation [27–30]. Both molecular dynamics methods and, to an even greater extent, quantum calculations require high computational capacity, which limits their application.

pacity, which limits their application. In our work, we have studied the subsequent reactivation steps using relatively simple and fast molecular docking methods that are available to a wide group of researchers. In addition to the reactivation of acetylcholinesterase, we have also explored this process In our work, we have studied the subsequent reactivation steps using relatively simple and fast molecular docking methods that are available to a wide group of researchers. In addition to the reactivation of acetylcholinesterase, we have also explored this process for butyrylcholinesterase. There is currently little research in the literature on this subject, although reactivation of BuChE appears to be beneficial from a pharmacological point of view. In addition, the presented studies consider an often overlooked state after reactivation when the OP fragment is bound to the oxime. The affinity of the oxime-OP complex to the

cholinesterase active site is important for the efficiency of the reactivation process. The present research broadens the knowledge about the cholinesterase reactivation process and provides methods and observations that can support future design of new AChE and BuChE reactivators.

#### **2. Materials and Methods**

The three-dimensional ligand structures (free oximes and oxime-OP complexes) were built with the Corina online tool (Molecular Networks GmbH, Nürnberg, Germany and Altamira, LLC, Westport, CT, USA). Free oximes were prepared in the anionic form. Sybyl-X 1.1 (Certara USA, Princeton, NJ, USA) was used for the assignation of formal and partial charges with Gasteiger-Marsili and the preparation of mol2 files with proper atom types.

All proteins used in the docking studies were prepared with Hermes 1.7 (Cambridge Crystallographic Data Centre (CCDC), Cambridge, UK) [31]. All histidine residues were protonated at Nε, and the hydrogen atoms were added. All complexes with modified OP conformations (5FPP-like conformation in the prereactivation state) or with OP fragments transferred from other structures (all OP-BuChE complexes) were previously optimized with the protein preparation tool from Maestro Suite (Schrödinger, New York, NY, USA) [32].

All docking experiments were performed with the Gold v5.3 program (CCDC, Cambridge, UK) [31]. A standard set of genetic algorithm parameters with a population size of 100 and a number of operations of 100,000 was applied. The binding site was always defined as all amino acid residues within 15 Å from the oxygen atom from the hydroxyl group of serine 203 (AChE) or 198 (BuChE). As a result, 10 ligand conformations were obtained and sorted according to the values of the following scoring functions: the Astex Statistical Potential (ASP), GoldScore, ChemScore and ChemPLP. All results were visualized by PyMOL 0.99rc2 [33].

Docking studies were divided into four parts.

In the first stage of this study, presenting a nonreactive conformation of oximes, ligands were docked to complexes of AChE (2Y2V and 3DL4) and BuChE (3DJY with OP transferred from AChE complexes) inhibited by sarin and tabun. No constraints were assigned to the time of docking. The optimal docking parameters identified during the tests for this stage are summarized in Supplementary Materials Table S1.

In the second stage, the previously used protein complexes were modified so that the conformations of the OP present corresponded to the sarin conformation observed in chain A of 5FPP. Sarin-AChE and sarin-BuChE complexes were obtained by transferring the sarin-bound Ser203 directly from the 5FPP complex and subsequently optimizing the structures to remove potential steric clashes. In the case of tabun-AChE, the OP was first covalently docked to the empty 3DL7 active site. Subsequently, the pose with an ethoxy substituent bent in a similar way to that of the sarin isopropoxy group from 5FPP was selected and optimized. The tabun-BuChE complex was built by replacing fragments that distinguished both organophosphates while maintaining the positions of the overlapping heavy atoms based on the 5FPP complex as a template. After transferring a serine-tabun fragment, the entire complex was optimized. To map the prereactivation position of the pyridinium oxime observed in the reference protein, substructure-based constraints were added with the oxime moiety from cocrystallization structure of asoxime with AChE (5FPP chain A). The most favorable docking parameters are summarized in Supplementary Materials Table S2.

In the third stage, all prepared oxime-OP complexes with the phosphorus atom in the trigonal bipyramidal geometry were covalently docked to the serine of the catalytic triad of each cholinesterase (AChE—PDB code: 1J06 and BuChE—PDB code: 1P0I), recreating the bond between the oxygen and phosphorus atoms. Selected docking parameters are summarized in Supplementary Materials Table S3.

During the last stage of this study, the OP-oxime (POX) complexes representing the postreactivation state (phosphorus atom in tetrahedral geometry) were docked to AChE (PDB code: 1J06) and BuChE (PDB code: 1P0I). Docking parameters allowing to obtain the most coherent results at this stage of research are summarized in Supplementary Materials Table S4.

#### **3. Results and Discussion**

For our research, we chose tabun and sarin as model organophosphorus compounds. The processes of reactivation of AChE and BuChE were analyzed for the following selected oxime compounds: 2-PAM, HI-6, obidoxime, K074, and K203. The ability of these compounds to reactivate AChE and BuChE blocked by sarin and tabun is presented in Table 1.



KD—phosphylated enzyme-oxime dissociation constant; kr2—second-order reactivation rate constant.

While modeling the reactivation process of cholinesterases blocked by organophosphorus compounds, we assumed the existence of several states in which the ligand could be found in the active site of the enzyme (Figure 3). Allgardsson et al. [19] proposed the presence of analogous states as part of the reactivation cycle. The crystal structure of AChE inhibited by sarin presented in their work shows a prereactivation state of HI-6 consistent with the mechanism of in-line attack of the reactivator on the phosphorus atom of an organophosphorus inhibitor [40]. As it is key for the reactivation process, we considered the balance between the formation of the presented complexes, which preferably should be shifted towards the stages immediately prior to reactivation. The energy differences between the prereactivation and transition complexes should be as small as possible, and the postreactivation OP-oxime complex should be poorly fitted to the active site. Moreover, to avoid additional enzyme blocking, the oxime should not show inhibitory effects. Free oximes were investigated as anions (monoanions in the case of obidoxime and K074) according to the reaction mechanism, indicating the direct participation of these forms in binding with OP. Poor dissociation of the oxime group was often indicated as the cause of the low activity of some reactivators in this analysis.

#### *3.1. Nonreactive Complex*

First, the nonreactive complex was analyzed. This complex is commonly observed in many crystals of AChE blocked by different OPs, e.g., in complexes stored under PDB codes: 5FPP (chain B), 5HFA, and 6CQV [20,41,42]. Our simulations also showed that the nonreactive complex is the conformation most readily achieved by ligands during docking to an active site blocked by OP. The analysis of this stage seems to be important both in the context of understanding the reasons for the behavior of the tested reactivators and to assess the preferences of individual ligands to adopt a prereactivation or nonreactive conformation. The information gained from this analysis will help guide further modification by shifting the balance between the two states in favor of the path leading to reactivation. To investigate the binding mode of the studied reactivators in a nonreactivating state, we docked them to the active site of AChE inhibited by sarin or tabun.

**Figure 3.** Multiple binding states during oxime reactivation of OP inhibited cholinesterase. **Figure 3.** Multiple binding states during oxime reactivation of OP inhibited cholinesterase.

*3.1. Nonreactive Complex*  First, the nonreactive complex was analyzed. This complex is commonly observed in many crystals of AChE blocked by different OPs, e.g., in complexes stored under PDB codes: 5FPP (chain B), 5HFA, and 6CQV [20,41,42]. Our simulations also showed that the nonreactive complex is the conformation most readily achieved by ligands during dock-After docking to the sarin-AChE complex, most poses obtained for 2-PAM revealed conspicuous cation-π and π-π stacking interactions with Trp286 or Tyr341 (Supplementary Materials Figure S1). In the case of the tabun-AChE complex, the trend was completely different. Most of the results obtained for 2-PAM revealed a dominating arrangement in which Trp86 was involved in interactions with the pyridinium ring through π-π stacking or cation-π interactions.

ing to an active site blocked by OP. The analysis of this stage seems to be important both in the context of understanding the reasons for the behavior of the tested reactivators and to assess the preferences of individual ligands to adopt a prereactivation or nonreactive conformation. The information gained from this analysis will help guide further modification by shifting the balance between the two states in favor of the path leading to reactivation. To investigate the binding mode of the studied reactivators in a nonreactivating state, we docked them to the active site of AChE inhibited by sarin or tabun. After docking to the sarin-AChE complex, most poses obtained for 2-PAM revealed Obidoxime revealed similar diversification in poses obtained for each AChE-OP complex (Supplementary Materials Figure S2). For AChE inhibited by sarin, the highest rated results showed the first pyridinium moiety interacting with Tyr124 (cation-π), Tyr337, and Tyr341 (CH-π). The other heterocyclic fragment was located in parallel with Trp286, allowing π-π stacking and cation-π interactions. AChE inhibited by tabun in most cases allowed obidoxime to interact with Trp86 by the pyridinium fragment in the parallel position while the other ring interacted with the PAS, most frequently with Trp286 and Tyr341.

conspicuous cation-π and π-π stacking interactions with Trp286 or Tyr341 (Supplementary Materials Figure S1). In the case of the tabun-AChE complex, the trend was completely different. Most of the results obtained for 2-PAM revealed a dominating arrangement in which Trp86 was involved in interactions with the pyridinium ring through π-π stacking or cation-π interactions. Obidoxime revealed similar diversification in poses obtained for each AChE-OP complex (Supplementary Materials Figure S2). For AChE inhibited by sarin, the highest rated results showed the first pyridinium moiety interacting with Tyr124 (cation-π), Tyr337, and Tyr341 (CH-π). The other heterocyclic fragment was located in parallel with Trp286, allowing π-π stacking and cation-π interactions. AChE inhibited by tabun in most Among the poses obtained for HI-6 during docking to sarin-AChE and tabun-AChE, we observed poses showing HI-6 with the amide group directed towards the catalytic site (Figure 4). A single pose for the sarin-AChE complex in a similar position to that from the 5FPP nonreactivated state was found. This dissimilarity resulted from the different conformations of Trp286, which, in our research, corresponded to the position observed in the unbound protein. Considering the above-mentioned single pose with an oxime fragment directed towards the interior of the active site of the enzyme, we observed an arrangement where the 1,2-disubstituted ring did not create any interactions, whereas oxygen belonging to the linker interacted with Tyr124 through hydrogen bond. The 1,4 disubstituted pyridinium ring interacted with Tyr341 (CH-π).

cases allowed obidoxime to interact with Trp86 by the pyridinium fragment in the parallel

Tyr341.

position while the other ring interacted with the PAS, most frequently with Trp286 and

we observed poses showing HI-6 with the amide group directed towards the catalytic site (Figure 4). A single pose for the sarin-AChE complex in a similar position to that from the 5FPP nonreactivated state was found. This dissimilarity resulted from the different conformations of Trp286, which, in our research, corresponded to the position observed in the unbound protein. Considering the above-mentioned single pose with an oxime fragment directed towards the interior of the active site of the enzyme, we observed an arrangement where the 1,2-disubstituted ring did not create any interactions, whereas oxy-

Among the poses obtained for HI-6 during docking to sarin-AChE and tabun-AChE,

*Biomolecules* **2021**, *11*, x FOR PEER REVIEW 7 of 21

substituted pyridinium ring interacted with Tyr341 (CH-π).

**Figure 4.** Docking results representing the nonreactivation conformation of the HI-6 (green) at the active site of AChE and BuChE (gray) blocked by sarin and tabun (purple). The key amino acids involved in the interactions with the reactivator are marked with yellow sticks. There is a clear difference in the binding of HI-6 in the active site of AChE and BuChE (the predominance of binding **Figure 4.** Docking results representing the nonreactivation conformation of the HI-6 (green) at the active site of AChE and BuChE (gray) blocked by sarin and tabun (purple). The key amino acids involved in the interactions with the reactivator are marked with yellow sticks. There is a clear difference in the binding of HI-6 in the active site of AChE and BuChE (the predominance of binding at the AChE PAS region and the BuChE anionic site).

at the AChE PAS region and the BuChE anionic site). The results for K074 and K203 coincided with those of the corresponding bispyridinium compounds obidoxime and HI-6. In the case of K074 (Supplementary Materials Figure S3), the best rated pose for AChE inhibited by sarin was arranged such that Tyr124, Tyr337, and Tyr341 were engaged in interactions with the compound. The pyridinium ring located closer to the catalytic site interacted with Tyr124 through cation-π binding, while Tyr337 created a CH-π bond. Tyr341 interacted with both pyridine rings through CH-π interactions. K074 docking to tabun-inhibited AChE revealed a binding mode close to that of obidoxime docking to tabun-inhibited AChE; the first pyridinium ring interacted with Trp86 (CH-π), while the other interacted with Tyr341 (CH-π). K203 (Supplementary Materials Figure S4) in sarin-AChE was arranged analogically to the arrangement of HI-6 docking to sarin- and tabun-inhibited AChE and was arranged in a manner in which a The results for K074 and K203 coincided with those of the corresponding bispyridinium compounds obidoxime and HI-6. In the case of K074 (Supplementary Materials Figure S3), the best rated pose for AChE inhibited by sarin was arranged such that Tyr124, Tyr337, and Tyr341 were engaged in interactions with the compound. The pyridinium ring located closer to the catalytic site interacted with Tyr124 through cation-π binding, while Tyr337 created a CH-π bond. Tyr341 interacted with both pyridine rings through CH-π interactions. K074 docking to tabun-inhibited AChE revealed a binding mode close to that of obidoxime docking to tabun-inhibited AChE; the first pyridinium ring interacted with Trp86 (CH-π), while the other interacted with Tyr341 (CH-π). K203 (Supplementary Materials Figure S4) in sarin-AChE was arranged analogically to the arrangement of HI-6 docking to sarin- and tabun-inhibited AChE and was arranged in a manner in which a 1,4-disubstituted ring with an amide group was directed towards the interior of the catalytic triad, while a 1,2-disubstituted pyridinium oxime was directed towards the exit of the narrow gorge. Only a single pose of K203 in complex with AChE inhibited by tabun revealed an arrangement where the aldoxime pyridinium remained at the catalytic site; for this pose, only a CH-π interaction with Tyr341 could be distinguished.

1,4-disubstituted ring with an amide group was directed towards the interior of the catalytic triad, while a 1,2-disubstituted pyridinium oxime was directed towards the exit of the narrow gorge. Only a single pose of K203 in complex with AChE inhibited by tabun revealed an arrangement where the aldoxime pyridinium remained at the catalytic site; for this pose, only a CH-π interaction with Tyr341 could be distinguished. There are currently a limited number of complexes of OP-BuChE with reactivators available in the PDB database. The structure closest to the described nonreactive state is the 4AXB crystal structure, in which 2-PAM is bound to the aged soman-BuChE complex. In docking studies, oximes easily adopt a nonreactive arrangement. This indicates a high preference for the tested compounds to occur in this conformation. It is stabilized mainly by the π-π and cation-π interactions of the pyridinium oxime moiety with Trp82. In the case of bispyridinium compounds, the other aromatic fragment created similar interactions with Tyr332.

Detailed analysis of the docking results for individual compounds revealed a very high similarity of the obtained poses both in the complex with sarin-BuChE and tabun-BuChE. In the case of 2-PAM (Supplementary Materials Figure S1), there was a slight difference in the arrangement of the pyridinium ring that directed the methyl substituent in a different direction. However, it had no significant effect on the formation of cation-π and π-π interactions with Trp82 or hydrogen bonds with the hydroxyl group of Tyr332.

Obidoxime (Supplementary Materials Figure S2) in both complexes created cation-π and π-π interactions with Trp82 and Tyr332 and hydrogen bonds between the ionized oxime moiety and Tyr128. At physiological pH, the reactivator, which has two oxime groups, was mostly present in the form in which one group was dissociated and the other was not. This allowed for the creation of an additional hydrogen bond between the nonionized oxime group and the Pro285 backbone in the sarin-BuChE complex, and, in the case of tabun-BuChE, a hydrogen bond with the carboxyl group of Glu197 was observed.

The compound HI-6 showed the greatest differences in its arrangement depending on the OP blocking the active sites of BuChE (Figure 4). In docking to sarin-BuChE, the most common and highest rated arrangement was a pose in which the pyridinium ring connected to the amide group formed cation-π and π-π interactions with Trp82. The amide group formed a hydrogen bond with Glu197. The pyridinium oxime fragment interacted with Tyr332 and was involved in the ionic bond with Asp70. The position of the oxygen atom from the linker indicates the possibility of creating hydrogen bonds with Tyr332. In the case of docking the compound HI-6 into tabun-BuChE, the results were significantly different. The pyridinium moiety with oxime substituents formed cation-π and π-π interactions with Trp82, as with most of the other examined reactivators. The oxime group was involved in creating a hydrogen bond with the hydroxyl group of Tyr332. In turn, the amide group formed a network of hydrogen bonds with the Leu286 and Ser287 main chains.

The results for compounds K074 and K203 (Supplementary Materials Figures S3 and S4) were very similar to those described for the structurally related obidoxime and HI-6. K074 at the active site of BuChE blocked by sarin created cation-π and π-π interactions with Trp82 through the ionized pyridinium oxime group. The other aromatic fragment created a CH-π interaction with Tyr332, and the oxime connected to it formed hydrogen bonds with Thr120 and the main chain of Ser287. In the tabun-BuChE complex, compound K074 was placed closer to the catalytic triad; therefore, K074 was the only tested reactivator that could form aromatic interactions with both Trp82 and His438. As in the case of obidoxime, we also observed the formation of hydrogen bonds between the ionized oxime group and Tyr128. The other aromatic ring formed cation-π and π-π interactions with Tyr332 and ionic interactions with Asp70. Compound K203 in the sarin-BuChE complex, similar to HI-6, interacted with Trp82 via a pyridinium ring with an amide substituent. Probably due to the longer linker, this fragment can be located closer to the catalytic triad so that the amide group forms hydrogen bonds with both Tyr128 and Glu197. The other pyridinium ring formed cation-π and π-π interactions with Tyr332 as well as ionic bonds with Asp70, and the oxime moiety was directed towards Thr120, forming a hydrogen bond with it.

The docking results indicate that in both AChE and BuChE, the aromatic rings of tryptophan residues greatly affect the arrangement of reactivators in the active site. In the case of AChE, the influence of Trp286 in the PAS seems to be predominant. In our research, we tried to limit the influence of amino acids that change their conformation because there is no precise research on the kinetics of these changes over time or their impact on reactivator binding. The obtained results indicate the need for further studies shedding more light on the contribution of Trp286. The change observed in the sarin-AChE complex described by Allgardsson et al. [19] significantly facilitates the positioning of the bispyridinium reactivators such that the oxime group is directed towards the interior of the active site of the enzyme. However, in BuChE, as the PAS is reduced and lacks this residue, interaction with Trp82 within the anionic site appears crucial. The interactions with the PAS that dominate in AChE enable the formation of both prereactivation and nonreactive

conformations. This indicates the ability to switch from one conformation to another. In turn, interactions with Trp82 in BuChE promote nonreactive conformations, which reduces the chance for effective reactivation of the enzyme.

#### *3.2. Prereactivation Complex*

The creation of a prereactivation complex is the first step in the process of reactivating serine blocked by an organophosphorus compound. This complex was experimentally confirmed in the 5FPP crystal structure (chain A) [19]. 5FPP represents HI-6 in the active site of sarin-inhibited AChE just before reactivation, in which the oxime group is directed towards the phosphorus atom, staying in close proximity. At the same time, this group remains in an axial relationship with the oxygen of Ser203. In addition to determining the position of the oxime at the moment of attack on the OP, a significant change in the position of the O-alkyl group is also observed. This change exposes the phosphorus atom, allowing the oxime anion to approach it.

To investigate this step of cholinesterase reactivation, modified crystal structures of AChE and BuChE inhibited by sarin or tabun were used. The modification concerned replacement of the original OP conformation with that known from 5FPP. This conformational change, which leads to exposure of the phosphorus atom to the attack of the oxime moiety, appears to be critical for obtaining the prereactivation state of the reactivator. During docking, we used constraints to obtain poses in which the oxime moiety is directed towards the OP phosphorus atom.

Visual inspection of the results obtained by docking to sarin-AChE showed two noteworthy clusters of poses. In the first cluster, the location of the oxime group resembled the one known from 5FPP, whereas the second cluster showed an aldoxime pyridinium ring shifted towards Trp86, allowing for π-π or cation-π interactions. There were also poses in which the aldoxime pyridinium ring was located at an intermediate position between those from the first and second clusters. The second cluster represented poses in which the aldoxime group did not turn towards the phosphorus atom of sarin, as the HI-6 in 5FPP also created bad contacts with the isopropoxy substituent of sarin. Accordingly, with the assumption that the oxime group remains in an axial relationship with Ser203 oxygen, we further described poses that belonged to the abovementioned first cluster of poses. In the case of 2-PAM (Supplementary Materials Figure S5), the best rated pose created cation-π and π-π stacking with Tyr341 and hydrophobic contacts with the acyl pocket. The 1,2-disubstituted ring of HI-6 (Figure 5) was slightly shifted towards Trp86 but remained between the phenol groups of Tyr124 and Tyr337, allowing cation-π stacking. In turn, the 1,4-disubstituted pyridinium moiety was parallel to Tyr341. K074 (Supplementary Materials Figure S7) interacted with Tyr341 and the acyl pocket as observed for 2-PAM. The pyridinium oxime fragment not directed towards the phosphorus atom was able to interact with Trp286 through CH-π stacking. Additionally, its oxime group interacted with the main chain amide moiety of Leu76. K203 showed π-π stacking interactions of the aldoxime pyridinium with Tyr341 and Phe338. The pyridinium ring with an amide moiety remained parallel at a close distance to Trp286. The amide group interacted with Tyr72 through hydrogen bonds. The obidoxime-binding mode (Supplementary Materials Figure S6) was the same as that for K203 (Supplementary Materials Figure S8), excluding interactions with Tyr72.

In the case of docking of reactivators to the tabun-AChE complex, for all scoring functions, the lowest rated oxime was 2-PAM (Supplementary Materials Figure S5), which also has the lowest reactivation ability. The scoring function results revealed accordance between their values and the activities of the reactivators. Oximes were arranged such that the oxime group was in axial relation with the Ser203 oxygen atom. To gain insight into pose clusters, three pose arrangements were distinguished. The first represents the group of oximes shifted towards the oxyanion hole, the second represents those shifted towards His447 and Phe338, and the third represents the intermediate position between the first and the second. Visual assessment of 2-PAM allowed us to observe interactions with Tyr124 and

region.

Tyr337 through cation-π and π-π stacking, while the oxime group was directed towards the phosphorus atom of tabun. The results obtained for HI-6 docking showed significant differentiation of poses (Figure 5). Based on the best rated poses, the key interaction seems to be cation-π interaction of the 1,2-disubstituted pyridinium ring with Tyr124. Additional contacts were made by the 1,2-disubstituted ring with Ser125 through cation-π interactions and a 1,4-disubstituted pyridinium ring with Tyr341 through cation-π and π-π interactions. Docking of K074 (Supplementary Materials Figure S7) revealed consistent poses that represented one binding mode. For example, the best rated pose for one pyridinium ring was placed at the catalytic site, creating cation-π interactions with Tyr124, while the other aromatic fragment created π–π stacking interactions with Trp286 and Tyr341 at the edge of the narrow gorge. The ionized oxime moiety was turned towards the amide group of the main chain of Arg296. Poses obtained in docking of K203 (Supplementary Materials Figure S8) showed high cohesion representing the same binding mode. The best rated pose showed the cation-π interactions of oxime containing a pyridinium ring with Tyr124, while the second amide ring created cation-π and π-π stacking interactions with Trp286. Furthermore, the amide group interacted with the main chain of Ser293 and Phe295 through hydrogen bonds. The best rated poses for obidoxime (Supplementary Materials Figure S6) were arranged as previously described for K074. Molecular modeling of the prereactivation state proved that the interaction of the pyridinium moiety with Tyr124 of AChE was crucial. *Biomolecules* **2021**, *11*, x FOR PEER REVIEW 10 of 21 S6) was the same as that for K203 (Supplementary Materials Figure S8), excluding interactions with Tyr72.

**Figure 5.** Docking results representing the prereactivation conformations of HI-6 (green) at the active sites of AChE and BuChE (gray) blocked by sarin and tabun (purple). The key amino acids involved in the interactions with the reactivator are marked with yellow sticks. The PAS region plays a huge role in ligand binding in both AChE and BuChE, despite the clear differences in the amino acids that build PAS of AChE and BuChE. Maintaining the BuChE prereactivation pose is dependent on the interaction with Tyr332 and Asp70. However, these interactions diminish during **Figure 5.** Docking results representing the prereactivation conformations of HI-6 (green) at the active sites of AChE and BuChE (gray) blocked by sarin and tabun (purple). The key amino acids involved in the interactions with the reactivator are marked with yellow sticks. The PAS region plays a huge role in ligand binding in both AChE and BuChE, despite the clear differences in the amino acids that build PAS of AChE and BuChE. Maintaining the BuChE prereactivation pose is dependent on the interaction with Tyr332 and Asp70. However, these interactions diminish during the formation of the transition state, in contrast to increasing interactions within the PAS AChE region.

the formation of the transition state, in contrast to increasing interactions within the PAS AChE After docking to BuChE, all the obtained prereactivation poses of the tested oximes created cation-π, π-π, or CH-π interactions with Tyr332. The ionic bond with Asp70 was also found in most results. Therefore, it can be assumed that both amino acids are crucial in

has the lowest reactivation ability. The scoring function results revealed accordance between their values and the activities of the reactivators. Oximes were arranged such that the oxime group was in axial relation with the Ser203 oxygen atom. To gain insight into pose clusters, three pose arrangements were distinguished. The first represents the group of oximes shifted towards the oxyanion hole, the second represents those shifted towards His447 and Phe338, and the third represents the intermediate position between the first and the second. Visual assessment of 2-PAM allowed us to observe interactions with Tyr124 and Tyr337 through cation-π and π-π stacking, while the oxime group was directed towards the phosphorus atom of tabun. The results obtained for HI-6 docking showed significant differentiation of poses (Figure 5). Based on the best rated poses, the key interaction seems to be cation-π interaction of the 1,2-disubstituted pyridinium ring with Tyr124. Additional contacts were made by the 1,2-disubstituted ring with Ser125 through cation-π interactions and a 1,4-disubstituted pyridinium ring with Tyr341 through cation-π and π-π interactions. Docking of K074 (Supplementary Materials Figure S7) revealed consistent poses that represented one binding mode. For example, the best rated pose for one pyridinium ring was placed at the catalytic site, creating cation-π interactions with Tyr124, while the other aromatic fragment created π–π stacking interactions with Trp286 and Tyr341 at the edge of the narrow gorge. The ionized oxime moiety was turned towards the amide group of the main chain of Arg296. Poses obtained in docking of K203 (Supplementary Materials Figure S8) showed high cohesion representing the

stabilizing the pyridinium oxime reactivators in the prereactivation position. However, the structure of the PAS in BuChE is not favorable for interactions with the studied compounds in this conformation. A significant difference between AChE and BuChE is the number of aromatic amino acids located near the entrance to the enzyme. When comparing the sequences of these proteins, we noticed that the aromatic residues building the PAS in AChE were replaced by short aliphatic or even polar amino acids in BuChE. Based on the 5FPP crystal structure, it is clear that for bispyridinium reactivators, the interactions with Tyr70, Tyr121, and Trp279, which are located in the PAS, are important for binding with AChE. In BuChE, these amino acids correspond to Asn68, Gln119, and Ala277, respectively. This is the reason for the worse fit of these reactivators to BuChE in the prereactivation conformation. The prereactivation position of 2-PAM (Supplementary Materials Figure S5) was generally the same in both the sarin-BuChE and tabun-BuChE complexes. The only specific observed interactions were cation-π and CH-π interactions with Tyr332.

The prereactivation arrangement of obidoxime (Supplementary Materials Figure S6) was also very similar for both complexes. The pyridinium ring with the attached ionized oxime moiety creates cation-π interactions with Tyr332 and ionic bonds with Asp70. In turn, the other pyridinium oxime fragment forms a hydrogen bond with the Leu273 main chain through the OH moiety. The position of this group can also be stabilized by hydrogen bonds with Asn68.

The differences in the positions of compound HI-6 were more significant than those for the previously described reactivators. Although the pyridinium oxime fragment was placed near Tyr332, creating the interactions previously described for the other compounds, the further aromatic fragment exhibited a more diverse arrangement. It formed hydrogen bonds with Asn289 (sarin-BuChE) or Glu276 (tabun-BuChE) through an amide moiety (Figure 5).

The results for K074 (Supplementary Materials Figure S7) and K203 (Supplementary Materials Figure S8) were similar to those for HI-6. K074 and K203 exhibited a number of interactions with Tyr332 (cation-π) and Asp70 (ionic). The position of the second heterocyclic fragment did not show such consistency. Among the results for K074 in the sarin-BuChE complex, it was observed that this fragment was directed towards the entrance to the enzyme, where it formed a hydrogen bond with Gln71. Concerning docking to tabun-BuChE, the pyridinium moiety was located deeper and formed hydrogen bonds with Gln119 and Asn289. K203, with a double bond in the linker, had reduced conformational flexibility, which resulted in a less favorable adjustment to BuChE near the entrance to the enzyme. In the case of the sarin-BuChE complex, there was a hydrogen bond between the amide group of K203 and Asn289. In the tabun-BuChE complex, the pyridinium ring substituted by the amide moiety of K203 did not create any specific interactions.

In summary, it is worth noting the similarities in the arrangement of reactivators in the cholinesterase active sites for this step. In both AChE and BuChE, the constraint-forced position of the oxime moiety promoted ligand conformation in which the aromatic fragment not involved in binding with the OP was directed towards Trp86/Trp82 (AChE/BuChE) or, for most of the results, towards the entrance to the enzyme active site. The conserved Asp74/Asp70, Tyr341/Tyr332, and Phe338/Phe329 residues significantly facilitated the pyridinium oxime fragment to adopt a prereactivation conformation. In AChE, there were additional aromatic residues in this area that considerably reduced space and thus oriented the reactivator directly to the OP-blocked serine. The larger space of the binding site in BuChE resulted in a shift of the pyridinium oxime moiety away from the position optimal for reactivation. This indicates the need to expand this fragment of the compound to compensate for the absence of aromatic residues that restrict the arrangement of the ligand. This modification could help the reactivator accommodate the prereactivation conformation at the BuChE active site. Another observation concerns the binding of the aromatic fragment within the PAS. As previously mentioned, the PAS in BuChE is significantly reduced compared to that in AChE. As a result, the specific interactions between the aromatic fragments of the ligands and the PAS in BuChE are considerably

weakened. The interactions are usually limited to hydrogen bonds formed by the amide or oxime moiety of the reactivator. The arrangement of this part of the ligand in BuChE often corresponds to that observed in AChE. However, in the case of the latter, the presence of additional aromatic residues, especially Trp286, that can adapt to the compounds stabilizes the beneficial conformation of the entire ligand. Therefore, it appears that to increase the efficiency of BuChE reactivators, the fragment of the compound not directly involved in the reactivation process needs to be optimized to better fit the reduced PAS in BuChE.

#### *3.3. Transitional State*

The application of covalent docking allowed us to reproduce the next step of the reactivation process, i.e., the formation of transitional serine-OP-oxime complex at the cholinesterase binding sites. This is the step with the highest energy, involving and it involves the spatial reorganization of the substitutes and the creation of a covalent bond between the phosphorus atom and the oxygen atom of the oxime group with the simultaneous breaking of the phosphorus-serine bond. The serine-OP fragment changes the geometry from tetrahedral to trigonal bipyramidal with both serine and oxime oxygen atoms in one axis [19,43]. For the purpose of this step, OP-oxime complexes with the appropriate geometry were prepared. Then, they were covalently docked to AChE and BuChE with the creation of a bond between the oxygen atom of the serine from the catalytic triad and the phosphorus atom of the OP. As in the previous steps, docking was conducted using the four scoring functions available in GOLD. The criteria for selecting the best conformation included both the coherence and the reliability of the mapped poses. As the correct poses were considered, those in which the free oxygen atom from the OP interacted with the oxyanion hole and the O-alkyl and alkyl substitutes (amino-alkyl in the case of tabun) were arranged similarly to those observed in known crystal complexes.

All the scoring functions enabled us to obtain the poses of the investigated oximes with attached sarin, in which the sarin binding mode was similar to that found in 5FPP. Regarding the pyridinium oxime fragment bound to sarin, the common feature of all conjugates regardless of scoring function was interaction with the phenol group of Tyr124. Additionally, this fragment interacted with Tyr341 and Phe338 through π-π stacking. The best scored pose for HI-6-sarin, in addition to interactions with Tyr341 and Phe338, was stabilized through cation-π and π-π stacking with Tyr337 (Figure 6). The best rated poses of K074, K203 and obidoxime revealed a binding mode similar to that of HI-6 (Supplementary Materials Figures S10 and S12). The pyridinium ring of these bispyridine compounds not bound to sarin mostly interacted with Trp286, creating cation-π and π-π stacking. Moreover, the amide groups of HI-6 and K203 were involved in hydrogen bonds with the phenol moiety of Tyr72.

The transitional complex for 2-PAM docked to tabun-AChE revealed the lowest values among the conjugates for all scoring functions (Supplementary Materials Figure S9). The values for the poses of the other oximes varied among all scoring functions and did not correlate with the experimental reactivation properties. Visual inspection of the best rated poses corresponding to the transitional state of 2-PAM-tabun showed a consistent binding mode of tabun that can be observed in original tabun-AChE complexes, for example, 3DL4. The phosphoryl oxygen interacted with the oxyanion hole through hydrogen bonds, whereas the dimethylamine group was directed towards the acyl pocket. Unlike the 3DL4 complex, the alkoxy substituent in this structure was bent in a manner similar to that of the sarin isopropoxy moiety found in 5FPP. The pyridinium ring neighbored Trp286, Tyr337, and Tyr124. The best rated pose showed cation-π interactions between the pyridinium moiety and Tyr337. Poses obtained for all scoring functions were coherent, but they differed from one another due to rotation changing affinities to neighboring residues. Covalent docking of HI-6-tabun to Ser203 had an effect similar to that of the earlier docking of 2- PAM-tabun regarding the arrangement of tabun (Figure 6). Most poses appeared in similar arrangements and represented the same binding mode. The 1,2-disubstituted pyridinium ring interacted with Trp86, Tyr337, and Phe338. The oxygen atom located at the linker

fragment created a hydrogen bond with Tyr124. The 1,4-disubstituted pyridinium ring directed towards the exit of the narrow gorge was arranged between Trp286 and Tyr341. The amide group remained without any interactions. *Biomolecules* **2021**, *11*, x FOR PEER REVIEW 13 of 21

**Figure 6.** Docking results reflecting the formation of the transitional state by the HI-6 (green) at the active sites of AChE and BuChE (gray) blocked by sarin and tabun (purple). The key amino acids involved in the interactions with the reactivator are marked with yellow sticks. The interactions of the reactivator with AChE PAS region, observed in the prereactivation state, were further enhanced during the formation of the transitional state. In the case of BuChE, the weakening of the interactions **Figure 6.** Docking results reflecting the formation of the transitional state by the HI-6 (green) at the active sites of AChE and BuChE (gray) blocked by sarin and tabun (purple). The key amino acids involved in the interactions with the reactivator are marked with yellow sticks. The interactions of the reactivator with AChE PAS region, observed in the prereactivation state, were further enhanced during the formation of the transitional state. In the case of BuChE, the weakening of the interactions with Tyr332 and Asp70 in the transitional state is not compensated by the PAS amino acids.

with Tyr332 and Asp70 in the transitional state is not compensated by the PAS amino acids. The transitional complex for 2-PAM docked to tabun-AChE revealed the lowest values among the conjugates for all scoring functions (Supplementary Materials Figure S9). The values for the poses of the other oximes varied among all scoring functions and did not correlate with the experimental reactivation properties. Visual inspection of the best rated poses corresponding to the transitional state of 2-PAM-tabun showed a consistent binding mode of tabun that can be observed in original tabun-AChE complexes, for example, 3DL4. The phosphoryl oxygen interacted with the oxyanion hole through hydrogen bonds, whereas the dimethylamine group was directed towards the acyl pocket. Unlike the 3DL4 complex, the alkoxy substituent in this structure was bent in a manner similar to that of the sarin isopropoxy moiety found in 5FPP. The pyridinium ring neighbored Trp286, Tyr337, and Tyr124. The best rated pose showed cation-π interactions between the pyridinium moiety and Tyr337. Poses obtained for all scoring functions were coherent, but they differed from one another due to rotation changing affinities to neighboring res-The results obtained in the covalent docking of the tabun complex with K074 and K203 showed more varied poses (Supplementary Materials Figures S11 and S12). There were notable differences regarding the tabun arrangement, which was divided into two clusters. The first cluster was consistent with the previously described position regarding the 2-PAM-tabun and HI-6-tabun complexes, while in the other cluster, we observed phosphoryl oxygen directed towards Glu202. Regarding the proximal 1,4-disubstituted pyridine ring of these oximes, it is notable that it interacts with amino acids such as Trp86 and Tyr337 and with the phenol group of Tyr124. The best rated pose of K074 showed catalytic site interactions of the pyridinium ring with Tyr337 through T-shaped stacking while the other aromatic fragment interacted with Tyr341 through parallel π–π and cation-π stacking. The oxime group not linked with tabun did not show any significant bonding. The best rated result for K203 was arranged such that one pyridinium ring linked with the phosphorus atom created T-shaped π-π stacking with Tyr337. The other pyridinium ring was parallel to Tyr341, allowing π-π and cation-π interactions. The amide group remained unbound. The investigation of obidoxime (Supplementary Materials Figure S10) revealed interesting binding modes. The phenol group of Tyr124 simultaneously interacted with both pyridinium rings. Moreover, the aromatic fragment not directly bound to the OP remained at the PAS of AChE and interacted with Tyr341.

idues. Covalent docking of HI-6-tabun to Ser203 had an effect similar to that of the earlier docking of 2-PAM-tabun regarding the arrangement of tabun (Figure 6). Most poses appeared in similar arrangements and represented the same binding mode. The 1,2-disub-In the case of BuChE complexes, detailed analysis of the results showed some evolution of interactions characteristic of the prereactivation complex. The interactions with Tyr332 were crucial for stabilizing the position of all tested compounds; however, they were limited to cation-π interactions, with the exception of 2-PAM, which additionally created

pyridinium ring directed towards the exit of the narrow gorge was arranged between

The results obtained in the covalent docking of the tabun complex with K074 and K203 showed more varied poses (Supplementary Materials Figures S11 and S12). There were notable differences regarding the tabun arrangement, which was divided into two clusters. The first cluster was consistent with the previously described position regarding the 2-PAM-tabun and HI-6-tabun complexes, while in the other cluster, we observed phosphoryl oxygen directed towards Glu202. Regarding the proximal 1,4-disubstituted pyridine ring of these oximes, it is notable that it interacts with amino acids such as Trp86 and Tyr337 and with the phenol group of Tyr124. The best rated pose of K074 showed

Trp286 and Tyr341. The amide group remained without any interactions.

CH-π interactions. The movement of compounds deeper into the BuChE binding site hindered the formation of ionic bonds with Asp70. At the same time, it enabled the creation of cation-π and π-π interactions with Phe329 for most compounds. 2-PAM (Supplementary Materials Figure S9) shows the strongest interactions with Phe329 and Tyr332 in both the sarin-BuChE and tabun-BuChE complexes. The presence of equivalent amino acids in the AChE active site (Phe338 and Tyr341) might be one of the reasons for the very similar reactivation activity of 2-PAM towards AChE and BuChE.

Another analyzed compound, obidoxime (Supplementary Materials Figure S10), in addition to the previously described interactions with Phe329 and Tyr332, also created significant π-π stacking between the pyridinium oxime moiety not involved in binding the OP and Trp82 as well as a hydrogen bond between the nonionized oxime moiety and Glu197. This is a significant difference compared to the prereactivation pose, in which this pyridinium fragment was directed towards the entrance to the BuChE active site.

Compound HI-6, although similar in structure to obidoxime, did not obtain similar poses (Figure 6). The transitional complex with sarin-BuChE showed hydrogen bonds between the amide group of HI-6 and the Gln67 side chain as well as the Asn68 main chain. In turn, in a tabun-BuChE complex, asoxime was directed towards the entrance to the active site, creating a hydrogen bond with Tyr332. It is likely that changing the position of the oxime from 4- to 3- hindered the binding of the aromatic ring near Trp82.

K074 (Supplementary Materials Figure S11) in both the sarin- and tabun-blocked BuChE complexes created a transitional state in which a pyridinium oxime fragment not involved in binding the OP was directed towards the entrance to the enzyme, creating a hydrogen bond with the Gly283 main chain. The stiffening of the linker by introducing a double bond and replacing one oxime moiety with an amide group made the compound K203 the only one of the tested reactivators that did not form any specific interactions with Phe329 in any of the complexes. The K203 transitional state (Supplementary Materials Figure S12) in sarin-BuChE indicated additional cation-π interactions between the amidebound pyridinium and Trp82. In the case of tabun-BuChE, this fragment was directed towards the entrance to the active site (in a pose analogous to that of compound K074), but no bonds that would stabilize this position were observed. It can be concluded that the ability of the pyridinium oxime moiety to create an optimal interaction with Phe329 and Tyr332 is key in the formation of the transitional complex by the studied reactivators. The substituents attached to the pyridinium ring not bound to OP should stabilize this position, as observed for obidoxime.

The results of covalent docking to the transitional state complexes of OP-AChE and OP-BuChE revealed further significant differences in the behavior of the tested reactivators. The highest rated poses for the 2-PAM transitional state with AChE and BuChE blocked by sarin or tabun presented a very similar position of the pyridinium fragment. This position created aromatic interactions with the Phe338/Phe329 and Tyr341/Tyr332 residues. In the case of AChE, an additional interaction with Tyr337 was observed. This is in accordance with the statement made in the previous step regarding the importance of interactions with these residues for the proper orientation of the pyridinium oxime fragment in a prereactivation state subsequently leading to a transitional state. Only the conformational changes of His447 and Phe338 present in the tabun-AChE complexes affected the position of the pyridinium ring, directing it towards Trp86. The results for the bispyridine compounds displayed high consistency of poses regarding docking to AChE and significant variation in the case of docking to BuChE. These discrepancies were again caused by the different structures of PAS in both enzymes. The aromatic residues of Tyr72 and Trp286 in AChE facilitated the arrangement of the aromatic fragment not directly bound to the OP within the PAS and thus positioned the entire ligand along the gorge of the enzyme. This binding mode was retained even for the 3DL4 complex, which exhibited the change in His447 and Phe338 conformations seen in tabun-AChE complexes. Moreover, the poses in which the ligand was placed over Phe338 clearly emphasize the limitations caused by this change in access to serine blocked by tabun. In the case of BuChE, the results can be divided

into two main conformations. The first conformation is similar to that observed for AChE; however, in this case, the interaction with Tyr332 is crucial. The other conformation shows the U-shaped bend of the compound so that the pyridinium fragment not connected to the OP can interact with Trp82. It appears that the type of *O*-alkyl substituent of the OP affects the conformation of the ligand to some extent. In the case of a larger *O*-isopropyl fragment of the sarin, a smaller number of poses showed interactions with Trp82. In turn, for the *O*-ethyl fragment of the tabun, this number was higher. Knowing the difference in the efficiency of the tested reactivators in relation to cholinesterases blocked by tabun and sarin, it seems that the strength of ligand interaction with Trp86/Trp82 inversely correlates with their reactivation potency. This observation is reflected in the participation of this residue in the formation of the nonreactive conformation.

#### *3.4. Postreactivation Complex*

After breaking the bond with serine, the organophosphorus fragment is transferred to the reactivator, forming a stable complex. From the point of view of the reactivation process, the OP-oxime complex (POX) should not have affinity to the binding site [19]. This would indicate a low potential of such a complex acting as an enzyme inhibitor and allow for its rapid dissociation from the binding site. For the purpose of our research, we prepared a number of oxime-OP complexes that we then docked to the AChE and BuChE active sites to examine their arrangement and possible inhibitory potential.

Among the poses obtained as a result of the docking of the sarin-oxime and tabunoxime complexes to AChE, two clusters can be distinguished. In the first cluster, the poses were arranged such that the OP moiety remained at the catalytic site. The other cluster revealed an arrangement in which the pyridinium moiety not connected with the organophosphorus moiety was located at the catalytic site, whereas the sarin or tabun group was directed towards the outside of the narrow gorge. The best scored poses of the sarin-2-PAM conjugate belonged to the first cluster. Visual inspection showed that the phosphoryl oxygen atom of sarin was directed towards the oxyanion hole, while the pyridinium ring of 2-PAM (Supplementary Materials Figure S13) interacted with the aromatic ring of Tyr124 through cation-π and π-π stacking. The best rated results of sarin-HI-6 obtained in docking to AChE could be divided into two clusters. In the first cluster, we observed that the phosphoryl oxygen atom of sarin was adjacent to the oxyanion hole. However, the OP moiety did not create any interactions. The oxime component of the complex created cation-π and π-π stacking with Tyr341 and Trp286 via the pyridinium ring (Figure 7). Other sarin-oxime conjugates revealed coherent poses that belonged to the second cluster (Supplementary Materials Figures S14–S16).

Visual cluster analysis allowed us to notice that the organophosphorus group of tabun-oxime complexes with more active reactivators (K203, K074, and obidoxime) was preferably directed towards the exit of the narrow gorge, whereas for less potent oximes, it occurred more frequently at the catalytic site (Figure 7). The best rated pose in docking of tabun-2-PAM to AChE (Supplementary Materials Figure S13) showed an arrangement in which the phosphoryl oxygen of tabun was directed towards the oxyanion hole while the dimethylamine substituent was located near the acyl pocket. Ethoxy substituents occupied space at the exit of the catalytic pocket. The pyridinium ring of the investigated complex neighbored Tyr341 and Trp286 and also created cation-π interactions with the phenol group of Tyr124. In the case of tabun-asoxime complex, the best rated pose was arranged such that the OP moiety remained at the catalytic pocket without any significant interaction. However, the 1,2-disubstituted ring perpendicularly neighbored Tyr124, allowing for some π-π interactions. The phenol group of Tyr124 was directed towards HI-6, creating a hydrogen bond with an oxygen atom from the linker. The 1,4-disubstituted pyridinium ring was parallel to Tyr341 at a distance of 4–5 Å. The amide group of HI-6 was close to the main chains of Phe295 and Ser293, allowing hydrogen bonding. In the case of the tabun-K074 complex (Supplementary Materials Figure S15), the OP fragment remained in the catalytic site without any favorable interaction, while the 1,4-disubstituted ring

attached to the OP fragment created cation-π stacking with the phenol group of Tyr124. Free aldoxime interacted with Gly342 through hydrogen bonds. Regarding the best rated pose for tabun-K203 (Supplementary Materials Figure S16), in which the OP component remained at the catalytic site, four interactions between pyridinium rings and AChE were noted. The aromatic fragment bound to the OP interacted with Tyr124 and Tyr341 via π-π stacking, while the phenol group of Tyr124 additionally interacted through cation-π stacking. The pyridinium ring with an amide moiety was parallel to Trp286 and allowed cation-π and π-π stacking. The OP component of the investigated conjugate did not create any significant interactions. The pyridinium fragment directly connected with the OP in the obidoxime-tabun complex (Supplementary Materials Figure S14) interacted with Tyr124 and Tyr341 by cation-π and π-π stacking. Furthermore, the close distance between the free aldoxime hydrogen atom and Ser293 main chain oxygen allowed hydrogen bonding. *Biomolecules* **2021**, *11*, x FOR PEER REVIEW 16 of 21

**Figure 7.** Docking results reflecting the behavior of the created sarin-HI-6 and tabun-HI-6 conglomerates (green) after reactivation of the active sites of AChE and BuChE (gray). The key amino acids involved in the interactions with the OP-reactivator conglomerate are marked with yellow sticks. The postreactivation complex binds more strongly at the active site of BuChE than AChE. **Figure 7.** Docking results reflecting the behavior of the created sarin-HI-6 and tabun-HI-6 conglomerates (green) after reactivation of the active sites of AChE and BuChE (gray). The key amino acids involved in the interactions with the OP-reactivator conglomerate are marked with yellow sticks. The postreactivation complex binds more strongly at the active site of BuChE than AChE.

Visual cluster analysis allowed us to notice that the organophosphorus group of tabun-oxime complexes with more active reactivators (K203, K074, and obidoxime) was preferably directed towards the exit of the narrow gorge, whereas for less potent oximes, it occurred more frequently at the catalytic site (Figure 7). The best rated pose in docking The results of docking to BuChE obtained for the various scoring functions were very consistent, especially for the highest rated poses. The key amino acids for binding OP-oxime complexes could be easily identified. These include Trp82, which participated in the binding of all tested postreactivation complexes, and Tyr332, which stabilized the binding mode of compounds with two pyridinium fragments (obidoxime, HI-6, K074, and K203) (Figure 7).

of tabun-2-PAM to AChE (Supplementary Materials Figure S13) showed an arrangement in which the phosphoryl oxygen of tabun was directed towards the oxyanion hole while the dimethylamine substituent was located near the acyl pocket. Ethoxy substituents occupied space at the exit of the catalytic pocket. The pyridinium ring of the investigated complex neighbored Tyr341 and Trp286 and also created cation-π interactions with the phenol group of Tyr124. In the case of tabun-asoxime complex, the best rated pose was A significant difference in the binding mode caused by the difference in the structure of the OP attached to the reactivator was observed for pralidoxime (Supplementary Materials Figure S13). For the sarin-pralidoxime complex, the pyridinium fragment was located among the aromatic rings of Trp82, Phe329, and Tyr332, creating hydrophobic and cation-π interactions. The oxygen atom of the sarin fragment additionally formed a hydrogen bond with Ser198 or His438. For the tabun-pralidoxime complex, in the most common pose, the pyridinium fragment was located between the Trp231 and Phe329 aromatic rings and also

lowing for some π-π interactions. The phenol group of Tyr124 was directed towards HI-6, creating a hydrogen bond with an oxygen atom from the linker. The 1,4-disubstituted pyridinium ring was parallel to Tyr341 at a distance of 4–5 Å. The amide group of HI-6 was close to the main chains of Phe295 and Ser293, allowing hydrogen bonding. In the case of the tabun-K074 complex (Supplementary Materials Figure S15), the OP fragment remained in the catalytic site without any favorable interaction, while the 1,4-disubstituted ring attached to the OP fragment created cation-π stacking with the phenol group of Tyr124. Free aldoxime interacted with Gly342 through hydrogen bonds. Regarding the best rated pose for tabun-K203 (Supplementary Materials Figure S16), in which the OP component remained at the catalytic site, four interactions between pyridinium rings and AChE were noted. The aromatic fragment bound to the OP interacted with Tyr124 and Tyr341 via π-π stacking, while the phenol group of Tyr124 additionally interacted through cation-π stacking. The pyridinium ring with an amide moiety was parallel to Trp286 and allowed cation-π and π-π stacking. The OP component of the investigated conjugate did not create any significant interactions. The pyridinium fragment directly connected with the OP in the obidoxime-tabun complex (Supplementary Materials Figure S14) interacted with Tyr124 and Tyr341 by cation-π and π-π stacking. Furthermore, the created hydrophobic and cation-π interactions. The oxygen atom of tabun was directed towards the oxyanion hole, but the ability to form hydrogen bonds in this position was limited.

Regarding the other reactivators, the differences caused by the type of attached OP were negligible. After binding to the OP, obidoxime (Supplementary Materials Figure S14) was arranged with a free pyridinium oxime fragment placed under the indole ring of Trp82, creating cation-π and π-π interactions. Nonreactive oxime group was in a convenient position to create hydrogen bonds with Tyr128 and Glu197. The sarin-bound pyridinium oxime fragment was placed parallel to the Tyr332 side chain, promoting cation-π and π-π interactions. This position was further stabilized by ionic bonds with Asp70. The sarin fragment facing the entrance to the active site created an additional hydrogen bond with Ser72.

The docking poses for the sarin-HI-6 and tabun-HI-6 complexes exhibited many interactions that were also observed in the previous obidoxime complexes (Figure 7). The amide-bound pyridinium ring formed cation-π and π-π interactions with Trp82. Additionally, the position of the ring allowed the creation of hydrogen bonds between the amide moiety and Thr120 and Tyr128. The main difference was the position of the sarin fragment, which was not directed towards the entrance to the active site of the enzyme but was located close to the catalytic triad. A hydrogen bond between the oxygen of the alkyloxy group and His438 was observed. A tabun fragment from an analogous complex was moved away from the catalytic triad; however, it was located deeper in the active site in this complex than in the tabun-obidoxime complex. The docking poses for the complexes of K074 and K203 (Supplementary Materials Figures S15 and S16) with both tabun and sarin were analogous to those observed for the docking of the obidoxime-OP complexes. The pyridinium ring substituted with an oxime (K074) or amide (K203) moiety was located parallel to the Trp82 side chain, forming cation-π and π-π interactions. This position was additionally stabilized by a hydrogen bond created between the oxime group and Tyr128. The second pyridinium ring was involved in cation-π and π-π interactions with Tyr332. Regardless of whether the compounds K074 and K203 formed complexes with sarin or tabun, the OP fragment was directed towards the entrance to the active site of BuChE, creating hydrogen bonds with the hydroxyl groups of the Ser72 and Gln71 main chains.

Comparing the obtained models representing the last stage of reactivation for AChE and BuChE, we observed that the differences in the trends of the arrangement of the compound depended on the structure of the reactivator and not on the type of OP or the particular cholinesterase. Both pralidoxime and HI-6 showed poor consistency of results with a slight predominance of poses in which the oxime-associated part of the postreactivation complex was directed inside the enzyme. On the other hand, reactivators such as obidoxime, K074, and K203 were arranged in a consistent manner in which the part associated with the OP was arranged within the outer part of the PAS. While the whole process of transition from the conformation in which the OP is located deep in the active site of a cholinesterase just after its bond with serine was broken to the conformation in which it is directed outwards requires further research on the molecular dynamics, it seems that stronger reactivators are characterized by this preference. It is also understood that the preference of the complex to leave the active site should positively correlate with the speed of this process and thus the restoration of the physiological functions of cholinesterase. The large number of poses in which we observed the OP-associated fragment facing outward tended to form an aromatic interaction of the pyridinium ring with Trp86/Trp82 so that the final pose was similar to the nonreactivating pose observed for the free reactivator docking.

#### **4. Conclusions**

In our study, we wanted to demonstrate the usefulness of a simple flexible molecular docking tool to model different steps in the reactivation process. We were pleased to find that by selecting various docking techniques, such as free flexible docking, flexible docking with distance constraints, and covalent docking, we obtained promising results

that allowed us to conduct a comprehensive analysis of OP–reactivator–cholinesterase interactions. Many groups involved in research on the reactivation of cholinesterases focus on an important stage in the reactivation process: the moment when the reactivator adopts a prereactivation pose [20,41,44]. Based on the reactivation cycle presented by Allgardsson et al. [20] in which the reactivator can adopt nonreactivation or prereactivation conformations at the active site of OP-blocked cholinesterase, we performed an interaction analysis for five known reactivators modeling both states at AChE and BuChE active sites blocked by sarin or tabun. For both AChE and BuChE, the nonreactive conformation was reached by the reactivator more easily, which may suggest an equilibrium shift between these positions in favor of a nonreactivating position. Analysis of interactions allowed us to identify the amino acids Phe338/Phe329 and Tyr341/Tyr332, which occur in both AChE and BuChE, as key to obtaining the prereactivation conformation. Many teams using the MD and QM/MM methods in their research have also reached similar conclusions [29,45]. Our results show that the amino acids composing the PAS of AChE are an important factor facilitating the prereactivation conformation of reactivators of this enzyme. Differences between AChE and BuChE in the structure of this fragment are the reason for the significantly lower reactivation efficacy of the tested reactivators for BuChE. The results we obtained are in line with the observations of other teams and the mapped analogous arrangements observed in the complexes available in PDB. Vyacheslav et al. [46], using the MD approach on the phosphorylated enzyme-reactivator complex, compared simulations starting from prereactivation and nonreactive positions and showed that the prereactivation (apical) position of the studied reactivators was rather stable and thus maintained over the MD trajectory, while in the nonreactive (side) position, the hydroxamic group left the active site shortly after the beginning of the simulation. Unfortunately, the small number of complexes representing some of the stages presented here, especially the prereactivation, transient and postreactivation states, makes it difficult to refer to the experimental data. This was especially troublesome when developing models for BuChE. Our study is one of the few that have conducted a comprehensive analysis of the course of reactivation of BuChE blocked by an OP. On the one hand, this limits our ability to refine the method based on the results published so far, but it is an important step towards explaining the reasons for the high resistance of OP-inhibited BuChE to known reactivators. The observed differences in the positions of the analyzed reactivators in the active sites of AChE and BuChE and the numerous similarities constitute a good starting point for proposing structural modifications to the known reactivators. Based on our observations, it seems that to improve the effectiveness of reactivators for BuChE, the key is to introduce groups that will stabilize the molecule in the PAS in a manner analogous to the interactions observed in AChE. This change could be achieved by introducing moieties that interact with the amino acids forming the acyl loop.

In summary, we believe that the protocol we have presented allows rapid and effective analysis of many important stages of the reactivation process. The presented approach may be part of advanced algorithms utilizing quantum mechanics/molecular mechanics (QM/MM) or molecular dynamics (MD) methods. After optimization of the ligand assessment process, the described method can provide an independent screening tool for predicting the potential efficacy of reactivation by new compounds.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2218-273 X/11/2/169/s1, Table S1: Detailed information about the docking parameters for Nonreactivating Conformation—Free Docking, Table S2: Detailed information about the docking parameters for Prereactivation Step—Docking with Substructure Based Constraints (Oxime Moiety from Chain A of 5FPP Complex, Table S3: Detailed information about the docking parameters for Transition State— Covalent Docking, Table S4: Detailed information about the docking parameters for Postreactivation State—Free Docking, Figures S1–S16: Visualization of docking results showing nonreactivating (S1–S4), prereactivation (S5–S8), transitional (S9–S12) and postreactivation (S13,S14) poses for 2-PAM, obidoxym, K074 and K203 reactivators.

**Author Contributions:** Conceptualization, J.J. and M.B.; methodology, J.J., J.K. and M.B.; investigation, J.J. and J.K.; writing—original draft preparation, J.J., J.K. and K.Ł.; writing—review and editing, J.J., J.K., K.Ł., B.M., K.M., Y.-S.J. and M.B.; supervision, M.B. and B.M.; funding acquisition, B.M., Y.-S.J. and K.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Polish National Centre for Research and Development, grant number V4-Korea 3/2018, by Czech Science Foundation, grant number GA18-01734S and by University of Hradec Kralove, Faculty of Science, grant number VT2019-2021.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This research was carried out as part of the international project Visegrad Group (V4)-Korea Joint Research Program on Chemistry and Chemical Engineering.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Impact of Sucrose as Osmolyte on Molecular Dynamics of Mouse Acetylcholinesterase**

**Sofya V. Lushchekina <sup>1</sup> , Gaetan Inidjel 2,3, Nicolas Martinez 2,3, Patrick Masson <sup>4</sup> , Marie Trovaslet-Leroy 5,**† **, Florian Nachon <sup>5</sup> , Michael Marek Koza <sup>2</sup> , Tilo Seydel <sup>2</sup> and Judith Peters 2,3,\***


Received: 6 November 2020; Accepted: 11 December 2020; Published: 12 December 2020 -

**Abstract:** The enzyme model, mouse acetylcholinesterase, which exhibits its active site at the bottom of a narrow gorge, was investigated in the presence of different concentrations of sucrose to shed light on the protein and water dynamics in cholinesterases. The study was conducted by incoherent neutron scattering, giving access to molecular dynamics within the time scale of sub-nano to nanoseconds, in comparison with molecular dynamics simulations. With increasing sucrose concentration, we found non-linear effects, e.g., first a decrease in the dynamics at 5 wt% followed by a gain at 10 wt% sucrose. Direct comparisons with simulations permitted us to understand the following findings: at 5 wt%, sugar molecules interact with the protein surface through water molecules and damp the motions to reduce the overall protein mobility, although the motions inside the gorge are enhanced due to water depletion. When going to 10 wt% of sucrose, some water molecules at the protein surface are replaced by sugar molecules. By penetrating the protein surface, they disrupt some of the intra-protein contacts, and induce new ones, creating new pathways for correlated motions, and therefore, increasing the dynamics. This exhaustive study allowed for an explanation of the detail interactions leading to the observed non-linear behavior.

**Keywords:** cholinesterase; osmotic stress; neutron scattering; molecular dynamics; MD simulations

#### **1. Introduction**

Investigations of the impact of osmotic stress on biological systems can play major role in the case of extreme environments, for instance, high salinity in oceans, the Dead Sea, and salty lakes. Osmotic stress may also be responsible for physio-pathological processes and the iatrogenic effects of injected drugs. High concentrations of osmolytes in protein solutions have important biotechnological and pharmaceutical applications. Moreover, osmotic pressure can serve as a tool to study particular parts of proteins, if access to pores or other confined regions is not possible otherwise. Such a strategy was chosen in the present case to shed light on the dynamics of the gorge of mouse acetylcholinesterase (mAChE). Acetylcholinesterase (AChE) is a key enzyme in the nervous system that terminates neurotransmission at central cholinergic synapses and at neuromuscular junctions by hydrolyzing the neurotransmitter's acetylcholine. The structures of cholinesterases (ChE) from different species are very close to each other with, for instance, identity scores between human AChE (hAChE) and mAChE of more than 82%. Due to the interest of cholinesterases as biopharmaceuticals for pre- and post-exposure treatments of organophosphorus poisoning [1], it is important to investigate the effect of sucrose as a protein structure protectant on the molecular dynamics of a model cholinesterase.

The first solved crystalline structure of a cholinesterase was that of *Torpedo californica* AChE (TcAChE) [2]. It showed that the catalytic active site is located at the bottom of a deep (20 Å) and narrow gorge (diameter of about 5 Å at the narrowest point, called "bottle neck"). This means that the substrate hydrolysis takes place in a mostly closed space virtually isolated from the bulk solvent. Along the passageway, the channel is so narrow that substrates would have no access to the active site if the enzyme were too rigid. However, Tai et al. [3] demonstrated. by means of molecular dynamics (MD) simulations, that half of the atoms in mAChE are participating in the so-called breathing mode that is periodically opening the gorge. The groups of McCammon [4–6] and Xu [7,8] undertook simulations of mAChE over 1 ns and 20 ns, respectively. The authors analyzed the dynamics and the possible opening of a so-called backdoor and other possible side channels in detail to shed light on movements in the gorge. Water molecules and co-solutes participate as well in the transport of the substrate to the enzyme catalytic center. Finally, water is the co-substrate of AChE-catalyzed hydrolysis of acetylcholine and other esters.

The experimental study of the motions in and out of the gorge is feasible by means of incoherent neutron scattering, which is sensitive to movements mainly of the nuclei of hydrogen (H) atoms, as the incoherent neutron scattering cross section of H is much higher than for any other nucleus. However, neutron scattering does not allow for distinguishing individual particles, which means that one observes the averaged motion of the H nuclei, and of molecular subgroups to which they are bound, assuming that they are representative for the dynamics inside the protein. In principle, incoherent neutron scattering would thus highlight parts of the sample by exchanging H against its isotope deuterium (D) and follow the dynamics of H belonging to the solution only, if the protein were deuterated. However, as so far it was not possible to deuterate ChEs, such an approach had to be excluded in this instance. Instead, the idea was to expose mAChE to sucrose unplugging all molecules from the gorge, to compare such sample to a native one in pure water and to extract information about motions of the solvent molecules only.

mAChE is a very fast enzyme, especially for a serine hydrolase, functioning at a rate approaching that of a diffusion-controlled reaction. The high speed of the enzyme is essential for the rapid functioning of cholinergic synapses, with a turnover of 103–10<sup>4</sup> s −1 . Therefore, we postulate that water dynamics within the gorge could be speeded up to facilitate transport of the substrate, eventually even beyond the normal diffusion rate of bulk water molecules. Such effects were observed, for instance, when water molecules were confined within a hydrophobic matrix [9]. The gorge could certainly be assumed to resemble such a confined environment.

Recently, we measured the dynamics of hAChE under high hydrostatic pressure, and our data suggests that this enzyme enters an intermediate folded state at 1.5 kbar, before being fully denatured at 3 kbar [10]. Osmotic pressure is believed to have different effects on reactions, as already discussed theoretically [11,12] and proven for butyrylcholinesterase (BChE) [13]. By using water-cosolvent mixtures, the enzyme hydration state can be modified, and thus, also the molecular dynamics and the transfer velocity of ligands and substrates. An osmolyte added to the solvent can either modify the viscosity only, as for instance glycerol (molar mass of 92 g/mol), or have a real osmotic effect, as for instance the bigger molecule of sucrose (molar mass of 342.3 g/mol). Sucrose is not able to penetrate the active site gorge of AChE and plays therefore the role of a semi permeable membrane by pumping the water molecules from inside the enzyme. Comparing the dynamics of AChE in the presence or absence of sucrose may hence provide information about the dynamics of water inside the gorge. Such information is extremely difficult to get by any other technique. Both effects, viscosity and osmotic pressure, can possibly have a great influence on the internal dynamics as well as the global

diffusion of the protein and induce changes in enzymatic activity. Experiments with enzymes in the presence of an osmolyte would allow us to complete our understanding of how the external conditions affect AChE dynamics and the water dynamics inside the gorge.

To better understand our results, we have undertaken molecular dynamics (MD) simulations of the same system at various sucrose concentrations. MD simulations give access to time scales very similar to those probed by incoherent neutron scattering, which allows us to undertake a direct comparison of dynamics and gain a better understanding of mechanisms at the atomic level.

#### **2. Materials and Methods**

#### *2.1. Sample Preparation*

Full cDNA of mAChE were inserted into pGS vector and expressed in Chinese hamster ovary cells (CHO-K1 cells). The transfected cells were selected and maintained in BioWhittaker® UltracultureTM medium (Lonza, Belgium) containing methionine sulfoximide (50 µM). The enzymes secreted into the culture medium, were purified using a procainamide affinity chromatography and were concentrated as previously described [14]. The gene produces a monomer, but the enzyme dimerizes spontaneously at the concentrations used here without inter-monomeric covalent binding. In its monomeric form, mAChE has a molecular mass of 66 kDa.

Activity measurements were carried out at 25 ◦C according to Ellman method [15] using 1 mM acetylthiocholine (ATC) as the substrate and 0.5 mM 5-5′ -dithio-bis (2-nitrobenzoic acid) (DTNB) in 0.1 M phosphate buffer pH 7.0.

The preparation of hAChE was previously described in [14]. mAChE was prepared exactly according to the same protocol. Briefly, mAChE was first dialyzed against 25 mM ammonium acetate dissolved in D2O (pD = 7.0) and then freeze dried (12 h) at 220 K under vacuum. This salt free protein powder, placed in an appropriated aluminum sample container, was dried for 12 h at atmospheric pressure over P2O<sup>5</sup> and weighed. This measured weight was the sample dry weight (h = 0 g D2O/g dry powder, denoted by g/g). For neutron experiments, the sample was hydrated by vapor exchange over D2O, at ambient temperature, in a desiccator, until a final water content of about 0.4 g/g was achieved. To verify that no loss of material had occurred, and that the hydration state was the same, the samples were weighed before and after the neutron scattering experiments and no losses were detected.

We prepared four samples of mAChE with 0%, 5%, 10%, and 15% of deuterated sucrose and hydrated them over D2O. We have chosen deuterated sucrose from Omicron Biochemicals, Inc., USA, to avoid a significant contribution from its H atoms to the incoherent neutron scattering signal. For the same reason, we used D2O for hydration instead of H2O. The proportions were chosen so that the weight of sucrose corresponded to the (*w*/*w*)% of mass of sucrose per mass of (sample + D2O) (see Table 1). We also prepared the buffer in D2O with 0 and 15 (*w*/*w*)% sucrose, corresponding to the extreme values and interpolating the intermediate values, as this is needed for data correction.


**Table 1.** Sample masses used for the experiments.

#### *2.2. Elastic Incoherent Neutron Scattering*

The protein was probed elastically and quasi-elastically by incoherent neutron scattering. We used three different spectrometers at the Institut Laue Langevin (ILL), France: the cold neutron time-of-flight spectrometer IN6 [16], the thermal backscattering spectrometer IN13 [17] and the high-resolution backscattering IN16 [18]. The spectrometers give access to molecular motions within the time windows of approximately 20 ps, 100 ps, and 1 ns, respectively. The experiments were conducted consecutively using the same samples.

The elastic incoherent neutron scattering (EINS) intensity is given within the Gaussian approximation by the dynamic structure factor at zero energy exchange:

$$S\_{el}(Q, \omega = 0 \pm \Delta E) \approx S\_0 \exp\left(-\frac{1}{3} \langle u^2 \rangle Q^2\right) \tag{1}$$

where <*u* <sup>2</sup>> is the average atomic mean square displacement (MSD). ω and *Q* are the exchanged energy and momentum in units of *h*¯, respectively, and ∆*E* is the half width half maximum (HWHM) of the instrumental energy resolution related to the time window through Heisenberg's uncertainty principle. The Gaussian approximation supposes that any atom can only undergo harmonic isotropic motions around its equilibrium position. The average MSD is related to the flexibility of the sample at a given temperature. For *Q* → 0, the approximation is strictly valid, and it holds up to <*u* <sup>2</sup>> *Q*<sup>2</sup> ≈ 1. The MSD can then be obtained for each temperature by the slope of the semi-logarithmic plot of the incoherent scattering function through: (, = 0 ± ) ≈ exp ൬− <sup>1</sup> 3 〈 ଶ 〉 ଶ൰ *ω*

$$\langle \mu^2 \rangle \approx -3 \frac{d \text{InS}\_{el}(Q, \omega = 0 \pm \Delta E)}{dQ^2} \tag{2}$$

The transmission values were measured on IN13 and were all above 0.9. Consequently, multiple scattering effects should not be taken into consideration for the data treatment. In order to obtain the scattered intensities of the sample, scattering from the empty sample holder and the buffer were subtracted. For the intermediate sucrose concentrations, we interpolated the curves of the buffer measured at 0 and 15 wt%. The data were normalized to a 2 mm thick vanadium slab, a totally incoherent scatterer, providing the relative detector efficiency and the instrumental resolution. Absorption correction was based on the correction formula of Paalman–Pings [19]. The complete data reduction was carried out using the LAMP software available at ILL [20]. → ≈ 〈 ଶ 〉 ≈ −3 (, = 0 ± ) <sup>ଶ</sup>

Quasi-elastic neutron scattering (QENS) is the part of scattering where small amounts of energy are exchanged between the neutrons and the sample, giving rise to a broadening of the elastic peak. The quasi-elastic structure factor

$$S(\mathcal{Q}, \omega) = e^{-\frac{1}{3} \langle u^2 \rangle \mathcal{Q}^2} \left[ A\_0(\mathcal{Q}) \delta(\omega) + \sum\_i A\_i(\mathcal{Q}) L(\Gamma\_i, \omega) \right] \tag{3}$$

contains the Debye–Waller factor *e* − 1 3 h*u* 2 i*Q*<sup>2</sup> , an elastic part, proportional to a delta-function in ω and representing the particles whose motions are not resolved within the instrumental setup. Furthermore, it includes a sum over Lorentzians, L, which describe different motional contributions. In theory, the sum can be extended up to an infinite number of contributions, but in practice a few Lorentzian curves are sufficient to describe the data reasonably well without over-fitting them. Here, we used two Lorentzian curves for the buffer and two Lorentzians for the sample, as more Lorentzians did not improve the fit quality. The amplitudes *A*0(*Q*) and *A<sup>i</sup>* (*Q*) depend on the momentum transfer *Q*. One can demonstrate [21] that the elastic incoherent structure factor (EISF) *A*0(*Q*) is given by the elastic intensity divided by the sum of elastic and quasi-elastic intensities at t → ∞. It vanishes except for *Q* = 0 if only long-range diffusive motions were present, and it is indicative of the geometry of the diffusional process otherwise. For the rotational diffusion of a particle, the EISF is unity at *Q* = 0 and falls to a minimum at a *Q* value which is inversely related to the radius of gyration of the rotating particle. (, ) = ି ଵ ଷ 〈௨ <sup>మ</sup>〉ொ మ ()() + ()( , ) ൩ ି భ య 〈௨ <sup>మ</sup>〉ொ మ ω → ∞

For data analysis, the structure factor must be convoluted with the instrumental energy resolution, which can be mimicked by a vanadium measurement:

$$S\_{\exp}(Q\_{\prime}\omega) = S(Q\_{\prime}\omega) \circledast S\_{\text{res}}(Q\_{\prime}\omega) \tag{4}$$

*Biomolecules* **2020**, *10*, 1664

When fitting the structure factor, it is possible to separate elastic and quasi-elastic parts. The Lorentzians are characterized by parameters of the motions, which can be extracted from their HWHM Γ*<sup>i</sup>* . Indeed, the form of the HWHM as function of *Q*<sup>2</sup> informs about the movements present in the sample. A pure translational diffusion in an unrestricted homogeneous medium corresponds to Brownian motion [22], whose characteristics are a linear behavior of the HWHM in *Q*<sup>2</sup> . The structure factor in this case reads:

$$S\_B(Q, \omega) = \frac{1}{\pi} \frac{D\_\Gamma Q^2}{\omega^2 + \left(D\_\Gamma Q^2\right)^2} \tag{5}$$

where *D<sup>T</sup>* represents the translational diffusion coefficient and Γ*<sup>B</sup>* = *D<sup>T</sup> Q*<sup>2</sup> its HWHM. Such a simple relation usually does not hold in the crowded medium of a macromolecule. The atoms diffuse, for instance, through a motion called translational jump-diffusion [23], where the HWHM does not increase linearly, but deviates from such behavior at higher *Q* values. The atoms perform oscillatory motions around their equilibrium positions for a time τ0. After that they diffuse for a time τ<sup>1</sup> by continuous diffusion and the process is then repeated. The structure factor can be described by a single Lorentzian with the HWHM Γ*JD* as:

$$
\Gamma\_{\rm JD} = \frac{D\_{\rm T}Q^2}{1 + D\_{\rm T}Q^2 \tau\_0} \tag{6}
$$

At small Q-values, the slope of Γ*JD* gives again the translational diffusion coefficient *D<sup>T</sup>* and the relation reduces to Brownian's law, i.e., Γ*JD* = Γ*B*. At large Q-values, the HWHM tends to the asymptotic value Γ<sup>∞</sup> = 1/τ0, τ<sup>0</sup> being called the average residence time.

When the motions are further restricted, they can result in rotational diffusion within a confined space [24]. The structure factor describing such movements reads:

$$S\_{\rm R}(\text{Q}, \omega) = A\_0(\text{QR})\delta(\omega) + \frac{1}{\pi} \sum\_{l=1}^{\infty} (2l+1) A\_l(\text{QR}) \frac{\Gamma\_l}{\omega^2 + \Gamma\_l^2} \tag{7}$$

where the diffusing atoms are limited to a sphere of radius *R*, Γ*<sup>l</sup>* = *D<sup>R</sup> l*(*l* + 1), *D<sup>R</sup>* being the rotational diffusion coefficient. Volino and Dianoux [24] found an analytical expression for the EISF:

$$A\_0(QR) = \left[\frac{\Im j\_1(QR)}{QR}\right]^2\tag{8}$$

where *j*1(*x*) = sin *x*/*x* <sup>2</sup> − cos *x*/*x* is the first-order spherical Bessel function. Bellissent-Funel et al. [25] expanded the model for the EISF by an immobile fraction *p*, where *p* denotes strongly bound protons. In that case, Equation (8) becomes:

$$A\_0(QR) = p + (1 - p) \left[ \frac{3j\_1(QR)}{QR} \right]^2 \tag{9}$$

Although the jump diffusion and diffusion in a confined space models were first developed for motions in fluids, such models were also successfully applied for crowded media [26].

#### *2.3. Molecular Dynamics Simulations*

The crystallographic structure of mAChE dimers PDB ID 1J06 [27] was used as a source of the protein coordinates. The missing loop 258–264 was restored using Modeller 9.14 [28], and the missing Lys496(A) and Arg493(B) side chains were reconstructed by means of a VMD *psfgen* module. Hydrogen atoms were added to construct the hydrogen bonding network by using Reduce software [29]. Water molecules recognized in the crystal structure were included in the model system and other non-protein molecules were removed. To meet experimental protein concentrations in the solution and to avoid interactions of the proteins with its copies in PBC, two dimers (four protein molecules) were

placed in solvation boxes. Clusters of sucrose molecules were built by means of VegaZZ software [30]. The total number of TIP3P water molecules and sucrose molecules added were adjusted to meet experimental concentrations. The number of molecules and sizes of the resulting systems are provided in Table S1 containing the electronic supplemental information (ESI).

MD simulations were performed with the NAMD 2.11 program [31] with CHARMM36 force field [32,33] at the Lomonosov-2 Moscow State University supercomputer center [34]. Periodic boundary conditions and PME electrostatics were applied. The systems were maintained at a constant temperature 300 K for productive runs and under pressure 1 atm (NPT ensemble) by using Langevin dynamics and Nosé-Hoover barostat. The initially prepared system underwent a 3000-step minimization.

As sucrose molecules were added to the system as a regular grid, systems with sucrose were subjected to a long pre-equilibration to optimize the solvents with 1 fs time-steps. For this, the coordinates of protein atoms and surrounding 4 Å water shells were fixed. In the first time, systems were subjected to 20 ns MD simulations at 400 K, and in the second time cooled down to 300 K during 20 ns. Then 10 ns unconstrained runs were performed for final equilibration. A radial distribution function was used to control the equilibration of the sucrose solution (see ESI, Figure S1). Then, 1 µs productive MD-runs at 300 K with 2 fs time-steps were performed. The analysis of obtained MD trajectories was performed using VMD [35] software and custom scripts. For calculations of the number of hydrogen bonds, a distance of 3.2 Å and a 40◦ cutoff angle were used. The calculations of the number of atom–atom contacts involved hydrogen atoms. Mean square displacement values were calculated from MD trajectories with different lag times τ as follows:

$$
\langle \mu^2(\tau) \rangle = \left[ RMSF(\tau) \right]^2 = \langle \left[ \overline{r}(t+\tau) - \overline{r}(t) \right]^2 \rangle \tag{10}
$$

where *RMSF* stands for root-mean square fluctuations.

Principle component analysis (PCA) and Cartesian PCA (cPCA) were performed using ProDy software [36], dihedral angle PCA (dPCA) was performed using Carma software [37]. cPCA was performed excluding the most mobile terminal residues 1–10 and 530–543. PCA methods allow us to determine and visualize the directions in space corresponding to the axes of highest mobility.

#### **3. Results**

The experiments were performed between June and August 2013 on the instruments IN6, IN13 and IN16 (DOI:10.5291/ILL-Data.8-04-707) and we repeated one measurement in August 2014 on IN13.

#### *3.1. Incoherent Neutron Scattering*

IN13 is a backscattering spectrometer giving access to averaged local motions within the time range of about 100 ps, e.g., it mainly enables the observation of sub-molecular motions of the amino acid side chains and small rotations at the protein's surface.

The four samples were scanned elastically within the temperature range from 270 to 310 K in steps of 5 K. After correction and normalization of the raw data, we first calculated the EINS of Equation (1) summed over all available *Q*-values. Such a quantity has smaller error bars due to the increased statistics and only depends on external parameters as temperature and the sucrose concentration. In Figure 1a, we present the summed intensities as function of wt% of sucrose at all measured temperatures. Higher temperature has the expected effect of lowering the intensities. In contrast, increased intensity depicts slower dynamics as more neutrons are scattered in a way to stay within the instrumental time window. All curves show a non-linear dependence on the sucrose concentration, with a maximum of intensity around 5 wt%. At 10 wt%, the intensities drop significantly, but present a pronounced dependence on temperature. In contrast, at 15 wt% all points seem to converge. We further used Equation (2) to extract the MSD of the same data and to plot them also as function of sucrose concentration (see Figure 1b). The MSD are mainly inversely proportional to the square of the summed intensities [38], so here we find a minimum at around 5 wt%. As the data have

lower statistics, the errors are bigger and the variations higher, especially at 15 wt% where the protein quantity and therefore the data quality was below the other curves.

**Figure 1.** Samples measured as function of sucrose wt% for temperatures ranging from 270 to 310 K. (**a**) Summed elastic intensities. The error bars were calculated but are so small that they are within the symbols. (**b**) mean square displacement (MSD) extracted through Equation (2).

Figure 2 presents the same MSD extracted from IN13 data now as function of the temperature for different concentrations of sucrose. All samples follow a linear increase in temperature except the sample at the highest sucrose concentration, which presents rather high fluctuations due to low data quality. As the behavior at 5 wt% sucrose was unexpected, we reopened the sample holder, dried, and rehydrated the sample and repeated the measurement. The result is shown in the orange dashed curve, close to the red curve, except at low temperature, where we saw probably the melting of heavy water around 277 K the second time.

**Figure 2.** MSD extracted through Equation (2) from data taken on IN13 as function of temperature. The black curve corresponds to the sample with no sucrose, the red and orange curves to the sample with 5 wt% sucrose, the cyan curve to 10 wt% sucrose, and the green curve to 15 wt% sucrose.

IN6 has a time window of about 20 ps and gives access to smaller and faster motions, which are for instance directly influenced by the presence of water molecules. Here, we measured EINS and QENS spectra. Globally, the MSD are lower compared to those obtained on IN13, what is explained by the smaller time window of IN6, but followed the same trends (see ESI, Figure S2) confirming further the non-linear behavior around 5 wt% sucrose.

QENS data were collected at 310 K, the physiological temperature of living organisms, and analyzed with the formalism described above (see Equations (3) and (4)). First, the buffer (D2O in presence or not of sucrose) curves were fitted with two Lorentzian functions and their HWHM extracted. Water molecules can simultaneously undergo translational and rotational motions and, indeed, we obtained this type of motions. Figure 3 shows the HWHM of the first Lorentzian curve together with fits according to the jump-diffusion model described in Equation (6). The systematic dip around 2.5 Å−<sup>2</sup> can be ascribed to de Gennes narrowing consequent to coherent scattering. The second Lorentzian was about 5–6 times larger, corresponding to faster motions, and corresponded to rotational motions as identified by the absence of a Q-dependence (data not shown). − –

 **Figure 3.** Half width half maximum (HWHM) Γ of the one Lorentzian curve used to fit the Quasi-elastic neutron scattering (QENS) data of D2O in presence or not of sucrose taken on IN6 as function of Q<sup>2</sup> . The black curve corresponds to no sucrose, the red curve to 5 wt% sucrose and the cyan curve to 10 wt% sucrose. The lines are fits using the jump-diffusion model Equation (6).

τ <sup>−</sup> All four samples present the typical Q-dependence of water at ambient temperature [39] and in the presence of co-solutes, corresponding to random-jump diffusion and rotations. In particular, our curves compare very well to those presented by Grimaldo et al. [40] for D2O measured also on IN6 at the ILL, Grenoble. The addition of sucrose had the expected effect: viscosity was increased and thus the translational diffusion coefficient D<sup>T</sup> decreased from 1.17 (5) to 0.89 (4), 0.73 (4) and 0.61 (3) 10 −5 cm<sup>2</sup> /s with the increasing sucrose concentration. The residence times τ<sup>0</sup> obtained for the four cases were 3.3 (3), 5.2 (5), 5.9 (6), and 5.5 (5) ps, respectively. In comparison, pure bulk H2O water has the characteristics *D*<sup>T</sup> = (2.5 ± 0.1) 10 −5 cm<sup>2</sup> /s; and τ<sup>0</sup> = 1.1 ps [25], which represent a higher and lower limit for these parameters. The self-diffusion in D2O is moreover slightly below that in H2O [41].

−

For the analysis of QENS data of the complete samples, we fixed then the HWHM of the two Lorentzians attributed to the buffer as well as the ratio of the absolute intensities between them, but let their absolute values vary freely together. Then, two Lorentzians were added to account for the sample. Figure 4 presents an example of such a fit of the experimental data with one elastic peak, four Lorentzian curves and a straight line to account for the background, which gives very good fit results.

− χ **Figure 4.** Logarithmic plot of an example of a fit of experimental QENS data of sample in the presence of 5 wt% sucrose for Q = 0.9 Å−<sup>1</sup> at 310 K. As described above, two Lorentzian curves were used for the buffer, two Lorentzian curves for the sample, a delta function to account for the elastic peak and a straight line for the background, resulting in a fit quality factor χ <sup>2</sup> = 1.38.

representative of the radius of the sphere delimiting the particles' motions. Both quantities are We determined first the EISF A<sup>0</sup> (Q) for the three samples with no, 5 and 10 wt% sucrose (see Figure 5). According to the diffusion-in-a-sphere model and Equation (9), the value of the EISF at highest Q is representative for the proportion p of particles, e.g., H nuclei and molecular subgroups to which they are bound, seen as immobile within the instrumental time resolution. R is representative of the radius of the sphere delimiting the particles' motions. Both quantities are progressively increasing with sucrose concentration (see Table 2), although the radius is mainly the same for all samples within error bars. Due to the increasing viscosity of the solution in the presence of sucrose, more particles are slowed down at higher sucrose concentration.

**Figure 5.** Elastic incoherent structure factor (EISF) of the samples with no sucrose (black curve), with 5 wt% sucrose (red curve) and with 10 wt% (cyan curve) sucrose as function of Q.


**Table 2.** Proportion p of particles seen as immobile and radius R of the sphere.

Further, we determined the HWHM of the Lorentzian curves of the samples (see Figures 6 and 7), plotted them as function of Q<sup>2</sup> and fitted them with one of the models.

 **Figure 6.** HWHM Γ of the narrower Lorentzian curve used to fit the QENS data of mAChE in the presence or not of sucrose taken on IN6 as function of Q<sup>2</sup> . The black curve corresponds to no sucrose, the red curve to 5 wt% sucrose and the cyan curve to 10 wt% sucrose. The lines are fits using the jump-diffusion model Equation (6). − − τ

−

 **Figure 7.** HWHM Γ of the broader Lorentzian curve used to fit the QENS data of mAChE in the presence or not of sucrose taken on IN6 as function of Q<sup>2</sup> . The black curve corresponds to no sucrose, the red curve to 5 wt% sucrose and the cyan curve to 10 wt% sucrose. The lines are fits with a constant value corresponding to rotational diffusion.

Figure 6 presents the HWHM Γ of the narrower Lorentzian function, corresponding to larger movements. We reiterate here that we are probing the internal motional diffusion of small molecular subgroups rather than global diffusion of the protein in solution, which is beyond the scope of the instrumental resolution of IN6. In this case, the tendencies are inverted with the curve at 10 wt% being clearly above the two others, which are more similar in absolute values. The Q-dependence of the HWHM still resembles jump-diffusional motions, but the slopes and curvatures differ significantly for the three samples and we notices again the non-linear behavior. The translational diffusion coefficient DT, corresponding to the slope at small *Q*-values, varies from (2.9 ± 0.4) 10−<sup>5</sup> cm<sup>2</sup> /s without sucrose to (2.1 ± 0.3) 10−<sup>5</sup> cm<sup>2</sup> /s at 5 wt% and to (2.7 ± 0.2) 10−<sup>5</sup> cm<sup>2</sup> /s at 10 wt% sucrose. The residence times τ<sup>0</sup> obtained for the three cases were now (5.5 ± 0.2), (3.8 ± 0.2) and (2.3 ± 0.1) ps, respectively.

Figure 7 shows the HWHM of the broader Lorentzian curve, used to fit the QENS data. The obtained values are tainted with bigger error bars and do not show a systematic dependence on Q. These smaller molecular motions can be attributed to rotational diffusion, which does not depend on Q. The HWHM Γ is related to the rotational diffusion coefficient Γ<sup>R</sup> = 2D<sup>R</sup> (see Equation (7)). Here, the effect of sucrose has an opposite effect: an increase of the rotational diffusional coefficient from 0 to 10 wt% yielding the values for D<sup>R</sup> of (0.24 ± 0.02) meV, (0.36 ± 0.02) meV, and (0.4 ± 0.01) meV, respectively.

Finally, IN16 is the instrument with the highest energy resolution which thus allows us to probe motions up to about 1 ns. Such a time window might provide us with further information on local motions, however, it also includes motions of subgroups and domains. We recorded EINS data from 275 to 312 K and plotted the summed intensities (see ESI, Figure S3a) and MSD (see ESI, Figure S3b). At these longer times, the differences are becoming even more striking, but again the hierarchy between the samples is the same.

#### *3.2. MD Simulations*

Mean square displacement values h*u* 2 i calculated for protein over MD trajectories with different time windows (τ) showed nonlinear patterns, as it was observed experimentally: a decrease in mobility at 5 wt% sucrose, compared to 0 wt%, and an increase at 10 wt% sucrose solution (Figure 8a). This effect becomes more pronounced with the increase of τ. The per-residue RMSF distribution is ambivalent, showing increased and decreased mobility of some areas in all systems (see ESI, Figure S4).

PCA of the MD trajectories allows to separate low frequency motions and to visualize differences in mAChE mobility (Figure 8b–d). cPCA shows that in the 0 wt% sucrose sample, the most mobile segments are peripheral loops, including loop 258–264, not resolved in X-ray, and Ω-loop. In the 5 wt% sucrose solution, these movements are dampened (amplitude decreased). However, Ω-loop movements are accompanied by acyl loop movements due to the binding of sucrose molecules at the gorge entrance. In the case of 10 wt% sucrose, movements of most of α-helixes and loops constituting the gorge walls (Ω-loop and acyl loop) have considerably increased their variance compared to movements in more diluted solutions, while β-sheets spanning through the protein are not affected. Other representations including network views of the cPCA results are shown in the ESI, Figure S5. dPCA free energy landscapes have distinctive local basins in case of 0 wt% and 10 wt% sucrose solutions, while for 5 wt% free energy landscape domains are more connected (ESI Figure S5).

the gorge walls (Ω

〈 2 〉 time windows (τ) showed nonlinear patterns, as it was observed experimental

wt% sucrose solution, these movements are dampened (amplitude decreased). However, Ω

wt% sucrose, movements of most of α

– ray, and Ω

This effect becomes more pronounced with the increase of τ. The per

–

movements in more diluted solutions, while β

**Figure 8.** (**a**) MSD vs. τ for different sucrose concentrations. The first cPCA modes with the greatest variance and for mAChE in 0 wt% (**b**), 5 wt% (**c**) and 10 wt% (**d**) sucrose solutions.

mAChE conformational states corresponding to the deepest basins seen for the 10 wt% solution reflect movements of the Ω-loop, induced by interactions with the sucrose molecule. Figure 9a shows the interactions of the Ω-loop with the adjacent part of the protein. They are few: hydrogen bonds between Asp74 and Tyr341, Tyr77, and Lys348 main chain, Trp86 and Tyr449; π-cation interaction between Phe80 and Lys348 and T-stacking interaction between Phe80 and Trp439. Panel B shows that all these interactions are disrupted during MD simulation of mAChE in 10 wt% sucrose solution. Yet, the Ω-loop is only partially detached through intermittent interactions, catalytic triad remains operative, and Trp86 forms cation-binding site. Furthermore, the snapshot shows a typical situation for MD simulations of ChEs in sucrose solutions: the sucrose molecule at the entrance of the gorge interacting with the Ω-loop and acyl loop (Figure 9b and ESI Figure S6). This results in the increased mobility of the acyl loop (Figure 8) and correlated motions between the Ω-loop and acyl loop residues (ESI Figure S6).

interacting with the Ω

the Ω

reflect movements of the Ω the interactions of the Ω

τ for different sucrose concentrations. The first cPCA modes with the great

and Lys348 main chain, Trp86 and Tyr449; π

Ω dashes show hydrogen bonds, orange dashes show π π and π after 1 μs MD trajectory **Figure 9.** Ω-loop constituting the wall of the mAChE active site gorge (purple) and its interaction with the adjacent part of the protein, yellow dashes show hydrogen bonds, orange dashes show π-π and π-cation interactions, sucrose molecules are shown with yellow surface representation. (**a**) X-ray structure, (**b**) snapshot after 1 µs MD trajectory of the 10 wt% solution system, connections between atoms shown in panel (**a**) as interactions are maintained. Acyl loop forming the other wall of the gorge is shown in green.

Ω Indeed, despite the observed increase in mobility of Ω-loop and acyl loop, there is no sign of partial unfolding, induced by the increase in sucrose concentration, or change in the gyration radius *R*g, or in the solvent accessible surface (SASA) (Table 3).


**Table 3.** Averaged parameters derived from 1 µs MD trajectories (±SD).

The number of water molecules inside the active site gorge is decreasing in the 5 wt% system due to osmotic stress [42]. The depletion of water molecules increases the mobility of the gorge residues (ESI, Figure S4b), but cPCA shows that these are high frequency motions and the catalytic triad remains operative. In the 10 wt% sample, this number is increased due to the instability induced by sucrose molecules on the Ω-loop and acyl loop. It opens new windows including a so-called backdoor [43] and side-door [44] for the influx of water molecules and reduces the accuracy of estimated water molecules inside the gorge (ESI Figure S6).

The averaged numbers of hydrogen bonds and atom–atom contacts have high variances; however, some tendencies could be noticed. Between the 0 wt% and 5 wt% sucrose systems, the number of protein–water interatomic contacts and hydrogen bonds slightly decreases. At the same time, the sum of the numbers of protein–water and protein–sucrose interatomic contacts and hydrogen bonds for the 5 wt% system is higher than the numbers of protein–water interatomic contacts and hydrogen bonds for the 0 wt% system. It can be interpreted as follows, when approaching the protein surface, the sucrose molecules rather interact with it through water molecules than to displace them. In addition, this could be the sign of a stiffer surrounding of the protein.

Comparing the 5 wt% and 10 wt% sucrose systems, the total number of hydrogen bonds and atom–atom interactions decreases as the percentage of sucrose is increased, and the replacement of water with sucrose molecules is more pronounced. Direct interactions of sucrose molecules with the protein surface does not lead to a significant decrease in protein intramolecular hydrogen bonding.

These differences between the nature of interactions between sucrose molecules and the protein in the 5 wt% and 10 wt% sucrose systems are illustrated by the example in Figure 10. In the 5 wt% sucrose system, the sucrose molecule interacts with mAChE residues through water molecules, on the contrary, in the 10 wt% sucrose system some interactions are through water molecules, others are direct. The salt bridge between Lys332 and Glu431 is broken, but the distant residue Asp396 is brought closer. This effect is similar to the one described above, the crosstalk between Ω-loop and acyl loop mediated by sucrose molecules. It can be seen as new contacts (correlated motions) in the network view of cPCA results (ESI Figure S5).

) snapshots after 1 μs MD trajectory showing the **Figure 10.** 5 wt% system (**a**) and 10 wt% system (**b**) snapshots after 1 µs MD trajectory showing the interaction of sucrose molecules with the same residues.

#### **4. Discussion**

– The role of viscosity on the sub-ns molecular dynamics of proteins has already been investigated in the past by means of incoherent neutron scattering and a clear impact was determined [45–47] of the presence of glycerol or glucose. As expected, co-solutes decreased the MSD of lysozyme, a small globular protein, in these cases above 150 K, but prevented the freezing even at low temperature. We were able to observe here similar behavior in the case of pure D2O buffer, which is more and more slowed down under the effect of sucrose (see Figure 3). The corresponding translational diffusion coefficients are all below the value known from the literature for bulk H2O [39]. In agreement with these findings, the residence times increased. The latter effect is indicative for a stronger confinement

form a layer at the protein's surface and protect the water layer, which permits a higher

wt% sucrose, the motion of the protein's outer part

effect [48]. The result found for the EISF of our samples points towards the same direction, as it increases at high Q values with increasing concentration of sucrose, meaning that less particles participate to the motions resolved by the instrument IN6 in the presence of the sugar although the radius of confinement is almost constant.

On the contrary, the results concerning the summed intensities and MSD extracted from three different spectrometers with increasing time window, IN6, IN13, and IN16, tell us a different story. As expected, the MSD increase with the time window, as more and more movements become visible at longer scales. The same findings hold for MD simulations as shown in Figure 8a. The MSD also increase with temperature, as the additional thermal energy permits larger amplitudes. However, the MSD equally increase with the concentration of sucrose, except for the sample at 5 wt% which presents a strong non-linear behavior. In particular, for IN6 we observe a cross-over of the curves corresponding to 0 and 5 wt% sucrose at around 290 K. Obviously, at this short time scale of 20 ps, the corresponding motions set in only at higher temperatures.

Several explanations are possible for such insights and only the MD simulations executed in parallel on the same kind of samples permitted to verify and better understand these effects:


QENS measurements show furthermore that not only molecular dynamics are reduced around 5 wt% of sucrose, but that the nature of the movements is also changed. Whereas the curves at 0 and 10 wt% of sucrose for the HWHM resemble to jump-diffusion motions, the curve at the intermediate sucrose concentration looks stiffer. This agrees with the fact that the modes at the protein's surface are dampened according to the simulations. The translational diffusion coefficients are indeed similar at 0 and 10 wt% but decreased in between. A re-structuration at the protein's surface and a strong coupling between the water layer and the sucrose matrix seems thus likely.

MD simulations are in good agreement with the experimental data and confirm the outstanding role of the interaction between water molecules, the sugar matrix and the protein surface. This allows us to understand the dynamics of the protein, e.g., the nonlinear behavior increasing over the time window. The analysis of MD trajectories, though partly speculative where values have high variances, allows us to suggest the following explanation: with the addition of sucrose at 5 wt%, the protein's nearest environment becomes more populated or stiff, reducing protein mobility as seen with neutron scattering. With an increase of the concentration, sucrose molecules replace more water molecules and have more direct interactions with the protein's surface. At 5 wt%, the water shell protects

the protein from sticky sucrose molecules, at 10 wt% this shell is depleted and this disturbs protein dynamics. Sucrose molecules penetrating the protein surface disrupt some of the intra-protein contacts, and induce new ones, creating new pathways for correlated motions. This effect is amplified by interactions between sucrose molecules, more frequent in more concentrated solutions.

It remains uncontested that the study of a protein in the presence of high sugar concentrations is not directly comparable to its native form. However, our observations shed light on the mechanisms of water fluxes into and out of the gorge, including the active site. Our investigations permitted us to conclude on complex dynamics of the Ω- and acyl loop. Other parts of the enzyme or polypeptide segments may also display peculiar dynamics. This opens new windows for understanding the extraordinary catalytic power of cholinesterases. In particular, the influx of water molecules across a backdoor could be addressed by MD simulations in the presence and absence of sucrose. The functionality of a backdoor as it was suggested [43] remains a puzzling issue. In the present case, simulations allowed to get information about the number of water molecules inside the gorge for the different samples. In contrast to what we expected (e.g., a continuous decrease of water molecules due to the osmotic effect and a concomitant shrinking of the gorge) maybe up to a collapse, we observed a much more subtle interplay between the different components and an increase in protein mobility at the highest sucrose concentration. The latter can be ascribed to the instability induced by sucrose molecules on Ω-loop and acyl loop, which finally results in conformational perturbation of the gorge.

#### **5. Conclusions**

In the present work, we intended to investigate the dynamics of the enzyme mAChE around its narrow gorge and of the water molecules, which move quickly into and out of the gorge. As the per-deuteration of mAChE is not possible so far, we opted for a study in D2O solution and in the presence of increasing concentration of deuterated sucrose. The sugar exerts osmotic pressure on the water molecules in the gorge permitting to study these movements in more detail. Surprisingly, we found experimentally a non-linear behavior in the dynamics of the samples with increasing sucrose concentration, inducing first a slowing down of the dynamics at 5 wt% and then again a raise. Moreover, the nature of the dynamics was also found different at the minimum of mobility. Extensive MD simulations allowed us to study in detail the interactions of the different components leading to such behavior. Sucrose molecules interact with the surface of the protein and the entrance of the gorge at a lower concentration through the water layer, damping the motions at the surface, but increasing them inside the gorge. When increasing the sucrose concentration more, the sucrose molecules replace some of the water molecules at the surface, permitting again to more water molecules to enter the gorge and opening simultaneously new pathways, among them the hypothesized backdoor to the gorge. These multiple mutual influences and interactions permit us to explain very satisfactorily the experimental findings.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2218-273X/10/12/1664/s1, Table S1: Sample proportions used for the simulations. Figure S1: MSD extracted through Equation (2) from data taken on IN6 as function of temperature. Figure S2: Summed intensities (a) and MSD (b) extracted through Equation (2) from data taken on IN16 as function of temperature. Figure S3: Protein surface shown for the snapshot in Figure 10B. Figure S4: cPCA results presented as network view and with color-coded backbone mobility.

**Author Contributions:** Conceptualization, P.M., M.T.-L., F.N. and J.P.; Data curation, N.M.; Formal analysis, S.V.L., G.I. and N.M.; Funding acquisition, P.M., M.T.-L. and F.N.; Investigation, M.T.-L., F.N., M.M.K., T.S. and J.P.; Methodology, S.V.L. and J.P.; Project administration, J.P.; Software, S.V.L.; Supervision, J.P.; Writing—original draft, S.V.L. and J.P.; Writing—review & editing, J.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** F.N and M.T. were funded by the French Ministry of Armed Forces (Direction générale de l'armement et Service de Santé des Armées) under grant BIOMDEF PDH-2-NRBC-3-C-301. P.M.'s contribution to the biochemical part was supported by the Russian Science Foundation grant # 20-14-00155.

**Acknowledgments:** The authors thank the ILL for attributing an internship financing to G.I. and for beamtime on the various instruments. Computer modeling was carried out using equipment from the shared research facilities of the HPC computing resources at Lomonosov Moscow State University [34]. We acknowledge the Joint Supercomputer Center of the Russian Academy of Sciences for provision of computational time.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
