QM/MM Description of Newly Selected Catalytic Bioscavengers Against Organophosphorus Compounds Revealed Reactivation Stimulus Mediated by Histidine Residue in the Acyl-Binding Loop

Butyrylcholinesterase (BChE) is considered as an efficient stoichiometric antidote against organophosphorus (OP) poisons. Recently we utilized combination of calculations and ultrahigh-throughput screening (uHTS) to select BChE variants capable of catalytic destruction of OP pesticide paraoxon. The purpose of this study was to elucidate the molecular mechanism underlying enzymatic hydrolysis of paraoxon by BChE variants using hybrid quantum mechanical/molecular mechanical (QM/MM) calculations. Detailed analysis of accomplished QM/MM runs revealed that histidine residues introduced into the acyl-binding loop are always located in close proximity with aspartate residue at position 70. Histidine residue acts as general base thus leading to attacking water molecule activation and subsequent SN2 inline hydrolysis resulting in BChE reactivation. This combination resembles canonical catalytic triad found in active centers of various proteases. Carboxyl group activates histidine residue by altering its pK_a, which in turn promotes the activation of water molecule in terms of its nucleophilicity. Observed re-protonation of catalytic serine residue at position 198 from histidine residue at position 438 recovers initial configuration of the enzyme’s active center, facilitating next catalytic cycle. We therefore suggest that utilization of uHTS platform in combination with deciphering of molecular mechanisms by QM/MM calculations may significantly improve our knowledge of enzyme function, propose new strategies for enzyme design and open new horizons in generation of catalytic bioscavengers against OP poisons.

Introduction

De novo protein design based on methods of structural bioinformatics is an actively developing field of modern life science. New biocatalysts created with this approach may be used for various purposes including therapeutic intervention (Vellard, 2003). Given natural limitations on the number of recombinant enzymes’ variants that can be tested experimentally, computational methods provide unique opportunities to narrow polypeptide sequence space to a set with reasonable size to be evaluated in search for novel properties, such as: fitting of new amino acid sequences into the existing structures (Dahiyat and Mayo, 1997; Socolich et al., 2005); creation of new proteins with previously unknown topology (Kuhlman et al., 2003; Koga et al., 2012); generation of enzymes with altered substrate specificity (Chen et al., 2009); design of enzymes with completely new catalytic activity (Zanghellini et al., 2006; Siegel et al., 2010); optimization of enzymes’ binding affinity (Lippow et al., 2007; Tinberg et al., 2013); creation of artificial proteins capable of selective paired interaction (Grigoryan et al., 2009); improvement of enzyme thermostability (Korkegian et al., 2005) or immunogenicity (Osipovitch et al., 2012; Salvat et al., 2014). Computational de novo design of enzymes catalyzing chemical reactions for which there are no natural analogous (Jiang et al., 2008; Rothlisberger et al., 2008; Siegel et al., 2010) may be regarded as breakthrough opening new era in enzymology. Usage of computational methods can be applied for generation of virtual libraries of enzyme’s active center (Smirnov et al., 2016) or as a starting point for the subsequent rational design of new biocatalysts with improved activity toward substrate of interest (Zheng et al., 2014).

Butyrylcholinesterase (BChE), a highly abundant serum enzyme, is a suicidal inactivator of organophosphorus toxins (OP). The exact function of BChE is unclear, nevertheless it is regarded as potent natural antidote against acetylcholine esterase poisons (Ilyushin et al., 2013; Lockridge, 2015; Terekhov et al., 2015). Since BChE is required to be administered in stoichiometric or even larger amounts in comparison with OP, application of BChE as an antidote toward OP poisoning may be restricted. It is evident that variants of BChE capable of catalytic destruction of OP are highly demanded (Aharoni et al., 2004; Valiyaveettil et al., 2011). BChE double mutant G117H/E197Q showed both reactivation and slowing down of aging process. However, example of its practical application for the treatment of OP poisoning clearly demonstrated that OP binding was also significantly decreased. This way the inhibition of endogenous BChE and AChE by OP will occur much earlier than interception of OP by bioscavenger, and therefore no protective effect would be achieved (Malisi et al., 2009). In other words, the development of effective catalytic antidotes is possible only if they either have highly efficient binding (similar to the wild-type BChE) or have extremely effective catalysis (k₂/K_M≈10⁷ M^-1min^-1). That was demonstrated for mutant forms of paraoxonase and phosphotriesterase (Jbilo et al., 1994; Mesulam et al., 2002), which had a protective effect against 2 × LD₅₀ cyclosarin (GF) and VX when the enzyme was administered at a dose of only 0.2 and 2 mg/kg, respectively. All these findings clearly show the need for further development of design approaches to create new catalytic bioscavengers. Recently we utilized combination of calculations and high-throughput screening technologies in order to create new BChE variants capable for catalytic degradation OP pesticide paraoxon (Terekhov et al., 2017). However, precise reason for this emergence of reactivation ability is unknown. The purpose of this study is to determine the molecular mechanism of paraoxonase activity mediated by newly described BChE variants using hybrid quantum mechanical/molecular mechanical (QM/MM) calculations.

Materials and Methods

Recombinant BChE Variants Expression and Purification

Genomic DNA from clones 14 and 15 was extracted as described previously (Lõoke et al., 2011) and used as a PCR template with primer pair a/b: (a) GCTTACTCGGAAGATGACATCATAATTGCAACAA (b) CTCCTTCGATCGACGCATGCAGAAAGCTCTGG. The product was cloned into the pFUSE-BChE-6xHis plasmid using NheI and NcoI restriction sites. These plasmids were used for transfection of HEK-293F cells using 293fectin Transfection Reagent (Thermo Fisher Scientific) according manufactures’ recommendations. Expression of BChE variants was carried out in FreeStyle 293-F Cells media (Thermo Fisher Scientific) in 250 ml culture flasks for 5 days. The cultural media was collected, concentrated using 30 kDa cut-off filter (Millipore) and further purified using TALON Metal Affinity Resin (Clontech) according manufactures’ recommendations.

Enzyme Activity Assay

Bimolecular inhibition constant (k₁/K_i) of BChE variants was determined according to the Kitz–Wilson method (Kitz and Wilson, 1962). The residual BChE activity was determined according to the Ellman’s method (Ellman et al., 1961). The activity was measured using Varioscan Flash Multimode Reader (Thermo Fisher scientific) at λ_abs = 405 nm. All data was derived from three technical replicates.

Cluster Analysis

To generate conformation states of BChE variants we employed Rosetta loopmodel (Stein and Kortemme, 2013) approach and generated 100 structures for each variant using PDB ID 1XLW as a template, where DEP was removed due to software limitations (Nachon et al., 2005). To reduce dimension of possible conformers we clustered conformations with AffinityPropagation (Frey and Dueck, 2007; Reshetnikov et al., 2018) according to the position of histidine(s) residue(s) relative to the catalytic serine residue: (i) histidine(s) residue(s) should be directed toward catalytic serine residue; (ii) amino acid residues in acyl-binding loop do not overlap with DEP (after structural alignment with 1XLW). Due to the tautomeric nature of histidine residue initial structures were generated with both possible protonated forms. DEP residue was copied from 1XLW into every analyzed structure.

QM/MM Calculations

The molecular mechanics (MM) approach was applied to the MM part of the system using parameters from the parm99sb-ildn force field with corrections (Lindorff-Larsen et al., 2010). The quantum mechanics (QM) subsystem was described utilizing DFTB approach (Grimme et al., 2010, 2011; Gaus et al., 2012, 2013). Side-chain atoms of catalytic triad (Ser¹⁹⁸, His⁴³⁸, and Glu³²⁵), Glu¹⁹⁷, DEP, newly introduced histidine residues from acyl-binding loop and neighboring water molecules in radius of 5 Å from DEP and histidine residues were included in the QM part of the system. Each simulation system was filled with water molecules according to the TIP3P model, and the total charge was neutralized with Na⁺ or Cl^- ions. Water and ions were equilibrated around the protein-DEP complexes by carrying out a 10-ps MD simulation with freezed position of protein-DEP conjugate, 50-ps with restrained position and 1-ns without restraints. For each cluster 10 independent replicas of preparation runs were carried and then analyzed to pick only productive conformations (with any water molecule in direct access to DEP, which was assumed as attacking) for subsequent reaction. Prepared systems were subjected to the QM/MM simulation with the GROMACS/DFTB package (Abraham et al., 2015; Kubar et al., 2015). The time step was set to 0.2 fs. Temperature coupling with Velocity Rescale (Bussi et al., 2007) scheme allowed observation of the behavior of systems at 300K. The total length of each initial ranking simulation was set to 10 ps. A metadynamics approach was used to overcome the activation barrier (Laio and Parrinello, 2002). In simulation run we used one composite collective variables (CV): distance Og(Ser¹⁹⁸)–P(DEP) minus distance P(DEP)–Ow(attacking water). A Gaussian history-dependent repulsive potential with a width of 0.006 nm and height of 2 kJ/mol was added every 40 steps. Additional upper wall restraints with kappa constant of 3000 were put on values of d(Og–P) and d(P–Ow) higher than 2.7 Å (minimal distance observed in QM/MM equilibration runs) to keep system within reaction conditions. To simulate specifically in-line hydrolysis a restraint with kappa constant of 500 was applied on angle Og–P–Ow being 180°. PLUMED plugin was used to apply metadynamics and restraints (Tribello et al., 2014). In total 320 QM/MM metadynamics runs were performed.

After initial ranking thorough energy profile scan was carried on selected systems. For this we utilized state-of-the-art TTMetaD enhanced sampling approach tailored to specifically address calculations of chemical reactions (Dama et al., 2014; Sun et al., 2016, 2017). To achieve better performance QM subsystem was redefined to exclude unnecessary water molecules that do not participate in proton transfer. Asp⁷⁰ was added since it proved to be crucial during initial ranking simulations. We used the same CV as stated before with potential width of 0.0045 nm and height of 1 kJ/mol. TTMetaD bias factor was set to be 42 kT, and initial guess for transition states (TS) coordinates was derived from ranking round to be -0.04 nm and 0.04 nm for first and second TS, respectively. Each profile was explored in parallel by 28 separate walkers to produce cumulative sampling time around 400 ps which corresponds to 2 million steps. A set of additional potentials was imposed to qualify for reaction mechanism devised from ranking round and could be accessed in Supplemental Information. Trajectories were analyzed with python scripts with the help of MDAnalysis (Michaud-Agrawal et al., 2011), and snapshots of trajectories with minimal energy barriers were created with PyMOL package.

Higher-level QM single-point energy calculations were performed using PBE0 (Perdew et al., 1996; Adamo and Barone, 1998). To further improve accuracy of calculations we utilized D3BJ set of empirical corrections (Grimme et al., 2010, 2011). To speedup calculation we limited basis set to def2/TZVP (Weigend and Ahlrichs, 2005) basis set with excluded f orbitals and enabled RIJCOSX approximation (Kossmann and Neese, 2010).

Results and Discussion

Selection of BChE Variants Capable of Catalytic Inactivation of Paraoxon Utilizing Ultrahigh-Throughput Screening

BChE variants were selected by previously developed ultrahigh-throughput screening technology of individual clones with different functional activity based on the usage of microfluidic droplets of the double water-oil-water emulsion (Terekhov et al., 2017). Schematically this method is described on Figure 1A. The advantage of applied methodology consists in generation of a monodisperse emulsion, therefore all individual droplets have the same dimensions. In this case, the fluorescence signal generated by the individual drops depends only on the concentration of the encapsulated fluorescent product, which in turn is formed as a result of the enzyme reaction, occurred in the drop. In order to provide appropriate quantification of the enzyme activity, it is necessary to achieve its effective and reproducible expression on the cell surface. Therefore, the optimization of the genetic construct encoded BChE, anchored on the surface of the yeast cells, in terms of the promoter, the leader peptide and the anchor sequence was carried out (Figures 1B,D).

FIGURE 1

FIGURE 1. (A) The scheme of the universal ultrahigh-throughput screening (uHTS) platform for biocatalysts selection. For the selected biocatalyst (1) the gene library is constructed (2). The library is further cloned into the specialized vectors (3) and transformed into the Pichia pastoris yeast cells, that present biocatalysts on the cell surface (4). The cells with displayed biocatalysts are screened by uHTS platform: (5) single cells are compartmentalized in the aqueous droplets of a double w/o/w emulsion with the fluorogenic substrate utilizing microfluidic technology. (6) Droplets with cells expressing active enzyme and therefore containing the fluorescent product are sorted by FACS. (B) Optimization of the plasmid vector for cloning of the BChE library for expression of the enzyme on the yeast cell surface. Scheme of the genetic constructs for expression of enzyme exposed on the cell surface (top). Comparative analysis of the surface-displayed BChE expressed under control of different promoters (left panel) and amount of BChE displayed on the cell surface promoted by different polypeptide anchors (right panel). (C) Extracellular release of the RFP, fused with different signal peptides. The genetic construct is presented on top. (D) Staining of yeast cells transformed by genetic construct (shown on top) by polyclonal anti-BChE antibodies (green fluorescence). Red fluorescence corresponds to the intracellular expression of the mCherry. (E) The scheme of the optimized genetic construct coding for the BChE variants used for uHTS (top). Analysis of enzymatic activity of surface-displayed WT BChE by Ellman’s assay.

The initial plasmid vector contained DNA coding for red fluorescent protein (RFP) mCherry fused with HA-tagged BChE via “self-processing” F2A peptide (Yang S. et al., 2008). This construct provides simultaneous production of a reporter and enzyme using a single mRNA transcript. Inducible aldehyde oxidase 1 (AOX1) and formaldehyde dehydrogenase 1 (FLD1) promoters or constitutive glyceraldehydes-3-phosphate dehydrogenase (GAP) promoter were used to determine the effect of the promoter on the level of production of the biocatalyst on the surface of the yeast cell (Figure 1B, left panel). Analysis of intracellular and membrane-allocated fluorescence revealed that the most efficient expression was achieved in the case of AOX1 promoter.

The effect of the anchor sequence on the expression of the recombinant biocatalyst on the surface of the yeast cell was studied using plasmids that differ only in the anchoring gene: either containing the sequence of the agglutinin surface antigen 1 (SAG1) subunit or the C-terminal fragment of the sperm-egg adhesion 1 (SED1) protein. After 72 h of analytical expression, the cells were stained with green fluorescent antibodies to HA and analyzed by cytofluorometry (Figure 1B, right panel). Despite in the case of SED1 anchor sequence the amount of the enzyme was much higher, the SAG1 provided the narrow distribution of enzyme on the surface, which is of the fundamental importance for the creation of cell-based libraries. The minimum possible distribution of the amount of marshaled enzyme is necessary for the subsequent correct selection of improved biocatalysts from libraries, the choice of biocatalysts from a wide distribution can lead to the selection of molecules without improved properties, but simply from better producing cells. Thus, for the creation and further selection of libraries by the yeast cell display method, it is more expedient to use SAG1. SED1 can be used to obtain clones-producers with the maximum level of protein production (in the medium or on the surface of the cell). Finally, the genetic construct used for the BChE yeast display contained methanol-inducible promoter AOX1, sequences of the red fluorescent reporter protein mCherry, F2A peptide, HSA leader to transport the expression product to the culture medium, the gene sequence of the enzyme to be studied, the HA epitope for immunofluorescent detection of the anchored enzyme, and the SAG1 sequence is linked serine glycine linker, which provides the attachment of the enzyme to the surface of the yeast cell (Figure 1C).

In order to determine the efficiency of the leader peptides providing the secretion of the protein of interest to the external medium, the α and α-short, HSA, KILL, AMY, PHO, and GLU signal sequences (please refer Supplementary Methods for details) were implemented in the signal reading frame with the nucleotide sequence coding for RFP. After 72 h of analytical expression of 96 clones of each leader peptide variant, the cells were removed by centrifugation and the fluorescent signal in the supernatants was analyzed (Figure 1D). Leader peptides α-short and AMY provide a high level of release of the RFP, whereas HSA demonstrated variation in the level of release of the RFP by different clones. Final genetic construct for library generation contained FLAG-tagged BChE gene with α-short leader peptide and C-terminal SAG1 sequence linked through serine glycine linker under control of AOX1 promoter. Transformation of yeast cells by this plasmid resulted in the appearance of a fully functional protein on the surface of the yeast cell. The analysis of the number of anchored molecules and calculation of the specific activity of the enzyme revealed that the developed yeast cellular display allows to obtain around 8000 molecules on a surface with specific activity that equals to the 15% of the activity of BChE in solution (Figure 1E).

The resulting vector was used to create a cell-based library of BChE variants. An acyl-binding loop ²⁸⁴TPLSV²⁸⁸ was used for random mutagenesis. Two variants of BChE exhibiting resistance to inactivation by paraoxon were selected (Terekhov et al., 2017). We further compared activity of the selected variants and 10 randomly picked clones from the initial library against parental substrate butyrylthiocholine and paraoxon (Table 1). Our data suggested that BChE clones containing sequences ²⁸⁴HTIHG²⁸⁸ (clone 14) and ²⁸⁴PSHSG²⁸⁸ (clone 15) hydrolyzed paraoxon, whereas clone 14 retained a significant level of specific activity toward butyrylthiocholine. Before only one known BChE substitution G117H was known to grant reactivation ability (Lockridge et al., 1997). Deciphering the precise molecular mechanism of newly found reactivating variants can provide valuable insight into this emerged organophosphate activity and help to devise new strategies for design of catalytic bioscavengers against OP poisons.

TABLE 1

TABLE 1. Kinetic and structural analysis of specific esterase and paraoxonase activity of selected BChE variants.

Computational Analysis of Molecular Mechanism of Paraoxon Hydrolysis by BChE Variants

Presence of proline residue in initial acyl-binding loop in principle should resulted in distinct and rigid configuration of this loop. Substitution of proline residue may lead to the increased loop flexibility, which in turn should be accounted during modeling process. In order to address this issue we generated 100 possible loop conformations for selected clones 14 and 15. Computational analysis of every possible loop configuration would require enormous computational resources, thus we were forced to decrease quantity of initial models. Therefore, we performed cluster analysis of generated configurations and finally obtained 11 clusters for both BChE variants (Figure 2). To further reduce number of initial models we withdraw configurations in which side chain of histidine residues do not fit into the reaction hemisphere. Despite both mutants are capable for reactivation, their activity is significantly decreased in comparison with common enzymes. That in turn means that frequency of catalytic events, especially reactivation with water, is very rare. Absence of relevant data on whether one or two water molecules and which one of two histidine residues in clone 14 are involved in the catalysis brought additional difficulties in the modeling. We therefore utilized a two-step approach: first we performed a “wide” round of multiple independent and relatively short calculation runs to get insight on principle reaction paths and rank all starting systems. Then only for the best of them we performed “deep” round of calculations aimed at accurate reconstruction of energy landscape. For both rounds we utilized metadynamics technique to simplify overcoming of energetic barriers since, on the one hand, it does not require a priori knowledge of reaction mechanism and, on the other hand, gives high-quality free energy estimate if mechanism is specified (Sun et al., 2017). Combination of metadynamics with hybrid QM/MM simulation recently had been shown as a suitable instrument for studying of slow enzymatic reactions involving organophosphorus compounds (Smirnov et al., 2016).

FIGURE 2

FIGURE 2. Positions of histidine residue selected after calculation and clusterization of structures corresponding to the clone 14 (A) and clone 15 (B).

Analysis of ranking round trajectories with the lowest energy barriers for both BChE variants revealed that histidine residue coordinate attacking water molecule therefore acting as a general base (Figures 3B,G). We showed that two histidine residues in acyl-binding loop of clone 14 participate in reaction of paraoxon hydrolysis in completely different ways. Histidine residue at position 284 always acts as a general base, whereas His²⁸⁷ never does this and instead coordinates water molecule (Figures 3A–E). Protonation state of histidine residue acting as general base evidently is crucial for the reaction. The most successful runs for clone 14 with ranking barrier of 32–33 kcal/mol were observed in cluster 11 containing His²⁸⁴ with proton on the δ-nitrogen atom (Table 2). Alternative clusters or protonation states yielded much higher energetic barriers or totally disrupted the reaction. In case of clone 15 with the 30 kcal/mol ranking barrier of the beneficial configuration originated from cluster 11 with proton on the ε-nitrogen atom (Figures 3F–J and Table 3). In both clones reaction ends with the proton transfer from His⁴³⁸ to Ser¹⁹⁸ thus completely restoring the enzyme’s initial state (Figures 3E,J). Detailed analysis of the accomplished runs revealed that histidine residues are always located adjacent to the aspartate residue at position 70. This combination resembles canonical pair of histidine and either aspartate or glutamate residues in the catalytic core of different proteases (Buller and Townsend, 2013). Carboxyl group functions as an activator of histidine residue by altering its pK_a and therefore promoting deprotonation of corresponding nucleophile (Figures 3B,G). Similar positioning of histidine residue and activating amino acid was recently proposed as a rational strategy for design of BChE variants capable of reactivation (Lushchekina et al., 2018).

FIGURE 3

FIGURE 3. Computer analysis of the catalytic mechanism of paraoxon hydrolysis by selected BChE variants. (A–E) Snapshots of the reactivation reaction path of paraoxon hydrolysis by clone 14. (F–J) Snapshots of the reactivation reaction path of paraoxon hydrolysis by clone 15.

TABLE 2

TABLE 2. Results of QM/MM ranking analysis for BChE clone 14.

TABLE 3

TABLE 3. Results of QM/MM ranking analysis for BChE clone 15.

With all the valuable insight gained from ranking calculation round we note that energy barrier estimates carry only qualitative and relative value which originates mainly from the choice of CV—all proton transfer events are not explicitly treated as a part it—together with insufficient sampling time. Therefore for two best systems—cluster 11 of clone 14 with both histidines as δ-protonated and cluster 11 of clone 15 with ε-protonated histidine—we further reconstructed energy profiles in the second calculation round with dramatically increased sampling time. Obtained profiles clearly show five distinct states: reagents state, first transition state, intermediate, second transition state, and products state (Figure 4). Estimated energy barriers—24.12 and 23.71 kcal/mol for clone 14 and 15, respectively—are in good agreement with values of 23.06 and 22.96 kcal/mol calculated from experimental data using Eyring-Polanyi equation (Figure 5).

FIGURE 4

FIGURE 4. Proposed generalized reaction mechanism of paraoxon hydrolysis by clones 14 and 15 based on trajectory snapshots. RS, reactant state; TS1, first transition state; INT, intermediate; TS2, second transition state; PS, product state.

FIGURE 5

FIGURE 5. Energy profile for clones 14 and 15. Metadynamics profile shown in solid line. Reported are energy values for reaction states with PBE0 calculations shown in parentheses.

Worth to mention that while DFTB approach was already used to study enzymatic reaction profiles and reproduced structural features well (Gruden et al., 2017), the overall performance in energy profile reconstruction based on single point energy evaluation for stationary points was reported to be below desirable level (Kulakova et al., 2015; Vasilevskaya et al., 2016). However, it was mentioned that the lack of accounting for configurational entropy may cause significant errors in the free energy surface for enzymatic reactions (Yang Y. et al., 2008; Hayashi et al., 2012). Therefore we utilized DFTB in the context of TTMetaD calculations, which as a valuable addition provides free energy profiles of the entire reaction path and not just stationary points. On the other hand, since DFTB geometries are regarded as reliable (Gruden et al., 2017), we decided to further validate our findings by performing energy calculations with higher-level hybrid functional PBE0 which has proven to be both accurate giving typical error less than 2 kcal/mol and robust for different type of systems (Luo et al., 2011). Energy profile reconstructed this way on points derived from metadynamics calculations produced similar values for energy estimate—23.06 and 22.83 kcal/mol for clone 14 and 15, respectively—thus also demonstrating good agreement with the experiment and supporting the choice of our methods.

Summarizing, active center of combinatorially selected BChE variants mediates canonical SN2 inline hydrolysis of covalent complex with paraoxon by incorporating additional nucleophile-base-acid triad as an outcome of introduced mutations. Observed reprotonation of catalytic Ser¹⁹⁸ by proton transfer from His⁴³⁸ recovers initial state of the enzyme, ready for the next catalytic turn. Data reported herein suggest that BChE with multiple substitutions in principle may seize paraoxonase activity and gives new ideas on design strategies to finally produce an efficient organophosphate enzyme. Synergy of the uHTS platform and in silico QM/MM explanation and prediction of amino acid residues’ catalytic function may open new horizons in generation of bioscavengers against OP poisons including but not limited to paraoxon molecule.

Author Contributions

AlZ, ArZ, and AnG performed the QM/MM calculations. YM and TB designed the genetic construct. ST designed and performed the screening and analyzed the kinetic experiments. AS and MA performed the molecular cloning. OK and SP performed the kinetic experiments. AB wrote the paper. AlG and IS contributed to conception and design of the study and wrote the paper.

Funding

This work was supported by RSF Grant 16-14-00191. Computer resources were provided by the Research Computing Center of Moscow State University. Supercomputers “Lomonosov” and “Lomonosov-2” were used for the calculations. Experiments were partially carried out using the equipment provided by the IBCH core facility (CKP IBCH, supported by Russian Ministry of Education and Science, grant RFMEFI62117X0018).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2018.00834/full#supplementary-material

References

Abraham, M. J., Murtola, T., Schulz, R., Páll, S., Smith, J. C., Hess, B., et al. (2015). GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1–2, 19–25. doi: 10.1016/j.softx.2015.06.001

CrossRef Full Text | Google Scholar

Adamo, C., and Barone, V. (1998). Toward chemical accuracy in the computation of NMR shieldings: the PBE0 model. Chem. Phys. Lett. 298, 113–119. doi: 10.1016/S0009-2614(98)01201-9

CrossRef Full Text | Google Scholar

Aharoni, A., Gaidukov, L., Yagur, S., Toker, L., Silman, I., and Tawfik, D. S. (2004). Directed evolution of mammalian paraoxonases PON1 and PON3 for bacterial expression and catalytic specialization. Proc. Natl. Acad. Sci. U.S.A. 101, 482–487. doi: 10.1073/pnas.2536901100

PubMed Abstract | CrossRef Full Text | Google Scholar

Buller, A. R., and Townsend, C. A. (2013). Intrinsic evolutionary constraints on protease structure, enzyme acylation, and the identity of the catalytic triad. Proc. Natl. Acad. Sci. U.S.A. 110, E653–E661. doi: 10.1073/pnas.1221050110

PubMed Abstract | CrossRef Full Text | Google Scholar

Bussi, G., Donadio, D., and Parrinello, M. (2007). Canonical sampling through velocity rescaling. J. Chem. Phys. 126:014101. doi: 10.1063/1.2408420

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C. Y., Georgiev, I., Anderson, A. C., and Donald, B. R. (2009). Computational structure-based redesign of enzyme activity. Proc. Natl. Acad. Sci. U.S.A. 106, 3764–3769. doi: 10.1073/pnas.0900266106

PubMed Abstract | CrossRef Full Text | Google Scholar

Dahiyat, B. I., and Mayo, S. L. (1997). De novo protein design: fully automated sequence selection. Science 278, 82–87. doi: 10.1126/science.278.5335.82

CrossRef Full Text | Google Scholar

Dama, J. F., Rotskoff, G., Parrinello, M., and Voth, G. A. (2014). Transition-tempered metadynamics: robust, convergent metadynamics via on-the-fly transition barrier estimation. J. Chem. Theory Comput. 10, 3626–3633. doi: 10.1021/ct500441q

PubMed Abstract | CrossRef Full Text | Google Scholar

Ellman, G. L., Courtney, K. D., Andres, V. Jr., and Featherstone, R. M. (1961). A new and rapid colorimetric determination of acetylcholinesterase activity. Biochem. Pharmacol. 7, 88–95. doi: 10.1016/0006-2952(61)90145-9

CrossRef Full Text | Google Scholar

Frey, B. J., and Dueck, D. (2007). Clustering by passing messages between data points. Science 315, 972–976. doi: 10.1126/science.1136800

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaus, M., Cui, Q., and Elstner, M. (2012). DFTB3: extension of the self-consistent-charge density-functional tight-binding method (SCC-DFTB). J. Chem. Theory Comput. 7, 931–948. doi: 10.1021/ct100684s

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaus, M., Goez, A., and Elstner, M. (2013). Parametrization and benchmark of DFTB3 for organic molecules. J. Chem. Theory Comput. 9, 338–354. doi: 10.1021/ct300849w

PubMed Abstract | CrossRef Full Text | Google Scholar

Grigoryan, G., Reinke, A. W., and Keating, A. E. (2009). Design of protein-interaction specificity gives selective bZIP-binding peptides. Nature 458, 859–864. doi: 10.1038/nature07885

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimme, S., Antony, J., Ehrlich, S., and Krieg, H. (2010). A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 132:154104. doi: 10.1063/1.3382344

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimme, S., Ehrlich, S., and Goerigk, L. (2011). Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 32, 1456–1465. doi: 10.1002/jcc.21759

PubMed Abstract | CrossRef Full Text | Google Scholar

Gruden, M., Andjeklović, L., Jissy, A. K., Stepanović, S., Zlatar, M., Cui, Q., et al. (2017). Benchmarking density functional tight binding models for barrier heights and reaction energetics of organic molecules. J. Comput. Chem. 38, 2171–2185. doi: 10.1002/jcc.24866

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayashi, S., Ueno, H., Shaikh, A. R., Umemura, M., Kamiya, M., Ito, Y., et al. (2012). Molecular mechanism of ATP hydrolysis in F1-ATPase revealed by molecular simulations and single-molecule observations. J. Am. Chem. Soc. 134, 8447–8454. doi: 10.1021/ja211027m

PubMed Abstract | CrossRef Full Text | Google Scholar

Ilyushin, D. G., Smirnov, I. V., Belogurov, AA Jr, Dyachenko, I. A., Zharmukhamedova, TIu, Novozhilova, T. I., et al. (2013). Chemical polysialylation of human recombinant butyrylcholinesterase delivers a long-acting bioscavenger for nerve agents in vivo. Proc. Natl. Acad. Sci. U.S.A. 110, 1243–1248. doi: 10.1073/pnas.1211118110

PubMed Abstract | CrossRef Full Text | Google Scholar

Jbilo, O., L’hermite, Y., Talesa, V., Toutant, J. P., and Chatonnet, A. (1994). Acetylcholinesterase and butyrylcholinesterase expression in adult rabbit tissues and during development. Eur. J. Biochem. 225, 115–124. doi: 10.1111/j.1432-1033.1994.00115.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, L., Althoff, E. A., Clemente, F. R., Doyle, L., Rothlisberger, D., Zanghellini, A., et al. (2008). De novo computational design of retro-aldol enzymes. Science 319, 1387–1391. doi: 10.1126/science.1152692

PubMed Abstract | CrossRef Full Text | Google Scholar

Kitz, R., and Wilson, I. B. (1962). Esters of methanesulfonic acid as irreversible inhibitors of acetylcholinesterase. J. Biol. Chem. 237, 3245–3249.

PubMed Abstract | Google Scholar

Koga, N., Tatsumi-Koga, R., Liu, G., Xiao, R., Acton, T. B., Montelione, G. T., et al. (2012). Principles for designing ideal protein structures. Nature 491, 222–227. doi: 10.1038/nature11600

PubMed Abstract | CrossRef Full Text | Google Scholar

Korkegian, A., Black, M. E., Baker, D., and Stoddard, B. L. (2005). Computational thermostabilization of an enzyme. Science 308, 857–860. doi: 10.1126/science.1107387

PubMed Abstract | CrossRef Full Text | Google Scholar

Kossmann, S., and Neese, F. (2010). Efficient structure optimization with second-order many-body perturbation theory: the RIJCOSX-MP2 method. J. Chem. Theory Comput. 6, 2325–2338. doi: 10.1021/ct100199k

PubMed Abstract | CrossRef Full Text | Google Scholar

Kubar, T., Welke, K., and Groenhof, G. (2015). New QM/MM implementation of the DFTB3 method in the gromacs package. J. Comput. Chem. 36, 1978–1989. doi: 10.1002/jcc.24029

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhlman, B., Dantas, G., Ireton, G. C., Varani, G., Stoddard, B. L., and Baker, D. (2003). Design of a novel globular protein fold with atomic-level accuracy. Science 302, 1364–1368. doi: 10.1126/science.1089427

PubMed Abstract | CrossRef Full Text | Google Scholar

Kulakova, A., Lushchekina, S., Grigorenko, B., and Nemukhin, A. (2015). Modeling reactivation of the phosphorylated human butyrylcholinesterase by QM(DFTB)/MM calculations. J. Theor. Comput. Chem. 14:1550051. doi: 10.1142/S0219633615500510

CrossRef Full Text | Google Scholar

Laio, A., and Parrinello, M. (2002). Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 99, 12562–12566. doi: 10.1073/pnas.202427399

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindorff-Larsen, K., Piana, S., Palmo, K., Maragakis, P., Klepeis, J. L., Dror, R. O., et al. (2010). Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins 78, 1950–1958. doi: 10.1002/prot.22711

PubMed Abstract | CrossRef Full Text | Google Scholar

Lippow, S. M., Wittrup, K. D., and Tidor, B. (2007). Computational design of antibody-affinity improvement beyond in vivo maturation. Nat. Biotechnol. 25, 1171–1176. doi: 10.1038/nbt1336

PubMed Abstract | CrossRef Full Text | Google Scholar

Lockridge, O. (2015). Review of human butyrylcholinesterase structure, function, genetic variants, history of use in the clinic, and potential therapeutic uses. Pharmacol. Ther. 148, 34–46. doi: 10.1016/j.pharmthera.2014.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Lockridge, O., Blong, R. M., Masson, P., Froment, M. T., Millard, C. B., and Broomfield, C. A. (1997). A single amino acid substitution, Gly117His, confers phosphotriesterase (organophosphorus acid anhydride hydrolase) activity on human butyrylcholinesterase. Biochemistry 36, 786–795. doi: 10.1021/bi961412g

PubMed Abstract | CrossRef Full Text | Google Scholar

Lõoke, M., Kristjuhan, K., and Kristjuhan, A. (2011). Extraction Of genomic DNA from yeasts for PCR-based applications. Biotechniques 50, 325–328. doi: 10.2144/000113672

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, S., Zhao, Y., and Truhlar, D. G. (2011). Validation of electronic structure methods for isomerization reactions of large organic molecules. Phys. Chem. Chem. Phys. 13, 13683–13689. doi: 10.1039/c1cp20834a

PubMed Abstract | CrossRef Full Text | Google Scholar

Lushchekina, S. V., Schopfer, L. M., Grigorenko, B. L., Nemukhin, A. V., Varfolomeev, S. D., Lockridge, O., et al. (2018). Optimization of cholinesterase-based catalytic bioscavengers against organophosphorus agents. Front. Pharmacol 9:211. doi: 10.3389/fphar.2018.00211

PubMed Abstract | CrossRef Full Text | Google Scholar

Malisi, C., Kohlbacher, O., and Hocker, B. (2009). Automated scaffold selection for enzyme design. Proteins 77, 74–83. doi: 10.1002/prot.22418

PubMed Abstract | CrossRef Full Text | Google Scholar

Mesulam, M. M., Guillozet, A., Shaw, P., Levey, A., Duysen, E. G., and Lockridge, O. (2002). Acetylcholinesterase knockouts establish central cholinergic pathways and can use butyrylcholinesterase to hydrolyze acetylcholine. Neuroscience 110, 627–639. doi: 10.1016/S0306-4522(01)00613-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Michaud-Agrawal, N., Denning, E. J., Woolf, T. B., and Beckstein, O. (2011). MDAnalysis: a toolkit for the analysis of molecular dynamics simulations. J. Comput. Chem. 32, 2319–2327. doi: 10.1002/jcc.21787

PubMed Abstract | CrossRef Full Text | Google Scholar

Nachon, F., Asojo, O. A., Borgstahl, G. E., Masson, P., and Lockridge, O. (2005). Role of water in aging of human butyrylcholinesterase inhibited by echothiophate: the crystal structure suggests two alternative mechanisms of aging. Biochemistry 44, 1154–1162. doi: 10.1021/bi048238d

PubMed Abstract | CrossRef Full Text | Google Scholar

Osipovitch, D. C., Parker, A. S., Makokha, C. D., Desrosiers, J., Kett, W. C., Moise, L., et al. (2012). Design and analysis of immune-evading enzymes for ADEPT therapy. Protein Eng. Des. Sel. 25, 613–623. doi: 10.1093/protein/gzs044

PubMed Abstract | CrossRef Full Text | Google Scholar

Perdew, J. P., Burke, K., and Ernzerhof, M. (1996). Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865–3868. doi: 10.1103/PhysRevLett.77.3865

PubMed Abstract | CrossRef Full Text | Google Scholar

Reshetnikov, R. V., Stolyarova, A. V., Zalevsky, A. O., Panteleev, D. Y., Pavlova, G. V., Klinov, D. V., et al. (2018). A coarse-grained model for DNA origami. Nucleic Acids Res. 46, 1102–1112. doi: 10.1093/nar/gkx1262

PubMed Abstract | CrossRef Full Text | Google Scholar

Rothlisberger, D., Khersonsky, O., Wollacott, A. M., Jiang, L., Dechancie, J., Betker, J., et al. (2008). Kemp elimination catalysts by computational enzyme design. Nature 453, 190–195. doi: 10.1038/nature06879

PubMed Abstract | CrossRef Full Text | Google Scholar

Salvat, R. S., Parker, A. S., Guilliams, A., Choi, Y., Bailey-Kellogg, C., and Griswold, K. E. (2014). Computationally driven deletion of broadly distributed T cell epitopes in a biotherapeutic candidate. Cell. Mol. Life Sci. 71, 4869–4880. doi: 10.1007/s00018-014-1652-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, J. B., Zanghellini, A., Lovick, H. M., Kiss, G., Lambert, A. R., St Clair, J. L., et al. (2010). Computational design of an enzyme catalyst for a stereoselective bimolecular Diels-Alder reaction. Science 329, 309–313. doi: 10.1126/science.1190239

PubMed Abstract | CrossRef Full Text | Google Scholar

Smirnov, I. V., Golovin, A. V., Chatziefthimiou, S. D., Stepanova, A. V., Peng, Y., Zolotareva, O. I., et al. (2016). Robotic QM/MM-driven maturation of antibody combining sites. Sci. Adv. 2, e1501695. doi: 10.1126/sciadv.1501695

PubMed Abstract | CrossRef Full Text | Google Scholar

Socolich, M., Lockless, S. W., Russ, W. P., Lee, H., Gardner, K. H., and Ranganathan, R. (2005). Evolutionary information for specifying a protein fold. Nature 437, 512–518. doi: 10.1038/nature03991

PubMed Abstract | CrossRef Full Text | Google Scholar

Stein, A., and Kortemme, T. (2013). Improvements to robotics-inspired conformational sampling in rosetta. PLoS One 8:e63090. doi: 10.1371/journal.pone.0063090

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, R., Dama, J. F., Tan, J. S., Rose, J. P., and Voth, G. A. (2016). Transition-tempered metadynamics is a promising tool for studying the permeation of drug-like molecules through membranes. J. Chem. Theory Comput. 12, 5157–5169. doi: 10.1021/acs.jctc.6b00206

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, R., Sode, O., Dama, J. F., and Voth, G. A. (2017). Simulating protein mediated hydrolysis of ATP and other nucleoside triphosphates by combining QM/MM molecular dynamics with advances in metadynamics. J. Chem. Theory Comput. 13, 2332–2341. doi: 10.1021/acs.jctc.7b00077

PubMed Abstract | CrossRef Full Text | Google Scholar

Terekhov, S., Smirnov, I., Bobik, T., Shamborant, O., Zenkova, M., Chernolovskaya, E., et al. (2015). A novel expression cassette delivers efficient production of exclusively tetrameric human butyrylcholinesterase with improved pharmacokinetics for protection against organophosphate poisoning. Biochimie 118, 51–59. doi: 10.1016/j.biochi.2015.07.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Terekhov, S. S., Smirnov, I. V., Stepanova, A. V., Bobik, T. V., Mokrushina, Y. A., Ponomarenko, N. A., et al. (2017). Microfluidic droplet platform for ultrahigh-throughput single-cell screening of biodiversity. Proc. Natl. Acad. Sci. U.S.A. 114, 2550–2555. doi: 10.1073/pnas.1621226114

PubMed Abstract | CrossRef Full Text | Google Scholar

Tinberg, C. E., Khare, S. D., Dou, J., Doyle, L., Nelson, J. W., Schena, A., et al. (2013). Computational design of ligand-binding proteins with high affinity and selectivity. Nature 501, 212–216. doi: 10.1038/nature12443

PubMed Abstract | CrossRef Full Text | Google Scholar

Tribello, G. A., Bonomi, M., Branduardi, D., Camilloni, C., and Bussi, G. (2014). PLUMED 2: new feathers for an old bird. Comput. Phys. Commun. 185, 604–613. doi: 10.1016/j.cpc.2013.09.018

CrossRef Full Text | Google Scholar

Valiyaveettil, M., Alamneh, Y., Rezk, P., Biggemann, L., Perkins, M. W., Sciuto, A. M., et al. (2011). Protective efficacy of catalytic bioscavenger, paraoxonase 1 against sarin and soman exposure in guinea pigs. Biochem. Pharmacol. 81, 800–809. doi: 10.1016/j.bcp.2010.12.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Vasilevskaya, T., Khrenova, M. G., Nemukhin, A. V., and Thiel, W. (2016). Reaction mechanism of matrix metalloproteinases with a catalytically active zinc ion studied by the QM(DFTB)/MM simulations. Mendeleev Commun. 26, 209–211. doi: 10.1016/j.mencom.2016.04.010

CrossRef Full Text | Google Scholar

Vellard, M. (2003). The enzyme as drug: application of enzymes as pharmaceuticals. Curr. Opin. Biotechnol. 14, 444–450. doi: 10.1016/S0958-1669(03)00092-2

CrossRef Full Text | Google Scholar

Weigend, F., and Ahlrichs, R. (2005). Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: design and assessment of accuracy. Phys. Chem. Chem. Phys. 7, 3297–3305. doi: 10.1039/b508541a

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S., Cohen, C. J., Peng, P. D., Zhao, Y., Cassard, L., Yu, Z., et al. (2008). Development of optimal bicistronic lentiviral vectors facilitates high-level TCR gene expression and robust tumor cell recognition. Gene Ther. 15, 1411–1423. doi: 10.1038/gt.2008.90

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Yu, H., and Cui, Q. (2008). Extensive conformational transitions are required to turn on ATP hydrolysis in myosin. J. Mol. Biol. 381, 1407–1420. doi: 10.1016/j.jmb.2008.06.071

PubMed Abstract | CrossRef Full Text | Google Scholar

Zanghellini, A., Jiang, L., Wollacott, A. M., Cheng, G., Meiler, J., Althoff, E. A., et al. (2006). New algorithms and an in silico benchmark for computational enzyme design. Protein Sci. 15, 2785–2794. doi: 10.1110/ps.062353106

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, F., Xue, L., Hou, S., Liu, J., Zhan, M., Yang, W., et al. (2014). A highly efficient cocaine-detoxifying enzyme obtained by computational design. Nat. Commun. 5:3457. doi: 10.1038/ncomms4457

PubMed Abstract | CrossRef Full Text | Google Scholar