Next Article in Journal
o-Carboranylalkoxy-1,3,5-Triazine Derivatives: Synthesis, Characterization, X-ray Structural Studies, and Biological Activity
Next Article in Special Issue
Design of a New α-1-C-Alkyl-DAB Derivative Acting as a Pharmacological Chaperone for β-Glucocerebrosidase Using Ligand Docking and Molecular Dynamics Simulation
Previous Article in Journal
Biofilm Inhibitory Abscisic Acid Derivatives from the Plant-Associated Dothideomycete Fungus, Roussoella sp.
Previous Article in Special Issue
Investigations of Structural Requirements for BRD4 Inhibitors through Ligand- and Structure-Based 3D QSAR Approaches
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Targeted Pesticides as Acetylcholinesterase Inhibitors: Comprehensive Cross-Organism Molecular Modelling Studies Performed to Anticipate the Pharmacology of Harmfulness to Humans In Vitro

by
Milan Mladenović
1,*,
Biljana B. Arsić
2,3,
Nevena Stanković
1,
Nezrina Mihović
1,
Rino Ragno
4,5,
Andrew Regan
6,
Jelena S. Milićević
7,
Tatjana M. Trtić-Petrović
7 and
Ružica Micić
8
1
Kragujevac Center for Computational Biochemistry, Faculty of Science, University of Kragujevac, Radoja Domanovića 12, P.O. Box 60, 34000 Kragujevac, Serbia
2
Department of Mathematics, Faculty of Sciences and Mathematics, University of Niš, Višegradska 33, 18000 Niš, Serbia
3
Division of Pharmacy and Optometry, University of Manchester, Oxford Road, Manchester M13 9PT, UK
4
Rome Center for Molecular Design, Department of Drug Chemistry and Technologies, Faculty of Pharmacy and Medicine, Sapienza Rome University, P.le A. Moro 5, 00185 Rome, Italy
5
Alchemical Dynamics srl, 00125 Rome, Italy
6
School of Chemistry, University of Manchester, Oxford Road, Manchester M13 9PL, UK
7
Vinča Institute of Nuclear Sciences, University of Belgrade, P.O. Box 522, 11001 Belgrade, Serbia
8
Faculty of Sciences and Mathematics, University of Priština, Lole Ribara 29, 38220 Kosovska Mitrovica, Serbia
*
Author to whom correspondence should be addressed.
Molecules 2018, 23(9), 2192; https://doi.org/10.3390/molecules23092192
Submission received: 6 August 2018 / Revised: 24 August 2018 / Accepted: 24 August 2018 / Published: 30 August 2018
(This article belongs to the Special Issue Application of Computational Methods in Drug Design)

Abstract

:
Commercially available pesticides were examined as Mus musculus and Homo sapiens acetylcholinesterase (mAChE and hAChE) inhibitors by means of ligand-based (LB) and structure-based (SB) in silico approaches. Initially, the crystal structures of simazine, monocrotophos, dimethoate, and acetamiprid were reproduced using various force fields. Subsequently, LB alignment rules were assessed and applied to determine the inter synaptic conformations of atrazine, propazine, carbofuran, carbaryl, tebufenozide, imidacloprid, diuron, monuron, and linuron. Afterwards, molecular docking and dynamics SB studies were performed on either mAChE or hAChE, to predict the listed pesticides’ binding modes. Calculated energies of global minima (Eglob_min) and free energies of binding (∆Gbinding) were correlated with the pesticides’ acute toxicities (i.e., the LD50 values) against mice, as well to generate the model that could predict the LD50s against humans. Although for most of the pesticides the low Eglob_min correlates with the high acute toxicity, it is the ∆Gbinding that conditions the LD50 values for all the evaluated pesticides. Derived pLD50 = f(∆Gbinding) mAChE model may predict the pLD50 against hAChE, too. The hAChE inhibition by atrazine, propazine, and simazine (the most toxic pesticides) was elucidated by SB quantum mechanics (QM) DFT mechanistic and concentration-dependent kinetic studies, enriching the knowledge for design of less toxic pesticides.

Graphical Abstract

1. Introduction

Pesticides [1] are either chemical or biological agents, widely used in agriculture [2], that deter, incapacitate, kill, or otherwise discourage pests. Their inappropriate administration leads to contaminated food while intentional release pollutes the environment, especially the soil, groundwater and natural waterways [3]. The uptake of contaminated food and water is potentially toxic to humans and other species. Pesticides mostly, exert toxicity by inhibiting the enzyme acetylcholinesterase (AChE) [4] which degrades acetylcholine (ACh), an essential neurotransmitter in the central nervous system (CNS) of insects, rodents, and humans [5]. AChE competitive inhibitors interrupt the physiology of autonomic ganglia, as well as, neuromuscular, parasympathetic and sympathetic effector junctions [6], controlled by ACh. However, little is known about the pharmacology of pesticides acting as AChE inhibitors. This report considers the pharmacology of common agricultural pesticides used in particular in the Western Balkans region [7,8]: atrazine, simazine, propazine, carbofuran, monocrotophos, dimethoate, carbaryl, tebufenozide, imidacloprid, acetamiprid, diuron, monuron, and linuron (Table 1). Despite their widespread usage, experimental crystal structures deposited at Cambridge Crystallographic Data Centre (CCDC) are available only for simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12], highlighting the poor availability of either ligand-based (CCDC deposition) or structure-based (Protein Data Bank (PDB) deposition) information about these compounds.
Triazines, such as atrazine, simazine, and propazine [13], are used as herbicides in crabgrass supression. Carbamates (carbofuran and carbaryl) are used against potato beetle; carbofuran was found extremely toxic for humans [14]. Organophosphate pesticides (monocrotophos and dimethoate) are used as insecticides [15]. Dimethoate is widely used to kill insects and as an insecticide for crops, vegetables, orchards, as well as for residential purposes [16]. Tebufenozide (a dialkylhydrazine pesticide) is used for plant protection [17]. Imidacloprid and acetamiprid, as neonicotinoids, are prescribed for the control of sucking insect pests (aphids, whiteflies, plant-hoppers, and thrips) and a number of coleopteran pests [17,18]. Diuron (a methylurea pesticide) is used to protect fruit, cotton, sugar cane and wheat [13]. Linuron and monuron (phenylurea pesticides) are plant protection agents [19]. Lethal doses for house and field mouse Mus musculus, for all of the above-mentioned pesticides, are known (Table 1) [20,21,22,23,24,25,26,27,28,29,30,31]. Lethal doses for Homo sapiens were calculated (Table 1) according to the recommendations of Reagan-Shaw et al. [32], using Equation (1):
H E D = a n i m a l   d o s e   a n i m a l   K m , h u m a n   K m
where HED is human equivalent dose in mg/kg of body weight, animal dose is used dosage for the specific model organism in mg/kg of body weight, animal Km factor is equal to 3, human Km factor is equal to 37.
The overexposure to pesticides results in the AChE inhibition at brain and nervous system nerve endings, as well as with other types of AChE found in the blood [33]; as a consequence, the ACh concentration increases its effects. Although, the AChE inhibition by pesticides listed in Table 1 is firmly established and the lethal doses are known, their behaviour in the inter-synaptic space and their interaction within the AChE active site, are still poorly investigated. Bearing in mind that different ACh conformations can trigger different receptors (nicotinics and muscarinics, respectively), conformational analysis (CA, ligand-based (LB) approach) could be considered as a computational tool to predict pesticides behaviour before their interaction with AChE [34]. Due to the lack of crystal structures of any pesticides co-crystallized with AChE, CA can be used to gather pesticide conformations that are, upon the interaction with AChE, going to be transformed into the bioactive ones [35].
For either mouse or humans as a model organism, AChE is a homodimer (Figure 1a), formed of two independent units, comprising 614 amino acids [36]. The AChE active site comprises the aromatic anionic domain, which accommodates the ACh positive quaternary amine, and the esterase catalytic domain, containing the catalytic triad Ser203, His447, and Glu334 with the purpose of hydrolysing the ACh to acetate and choline. The ester hydrolysis reaction leads to the formation of an acetylated-enzyme. Then, the acetyl-enzyme undergoes nucleophilic attack by a water molecule, assisted by the His447, freeing the acetic acid, likely as an acetate ion, and regenerating the free enzyme (Figure 1b) [33].
Given that the pharmacology of the majority of widely used pesticides (Table 1) is still poorly understood, the focus of this study was to provide the LB or SB cheminformatics-based guidelines for several important tasks, and unresolved questions: (1) how to define the pesticides’ structures in aqueous solution prior to interaction with AChE? (2) how to define the pesticides’ bioactive conformations upon the interaction with AChE? Other important assignments of the LB or SB protocols applied herein were to provide answers to several important questions: (1) is there any possibility to anticipate the pesticides acute toxicity from pre-bound conformations, before the actual interaction with AChE? (2) can the pesticides’ acute toxicity be quantified from the bound conformations, upon the actual interaction with AChE? To the best of our knowledge, there are no previous reports on either LB or SB efforts to predict the bioactive conformations of investigated pesticides (Table 1) for both mouse and human AChEs, as well as to enumerate their acute toxicities by means of free energy of formation (Eglob_min) in the aqueous solution and/or free energy of binding (∆Gbinding) to AChE. All the cheminformatics findings were externally validated on the series of other commonly used pesticides (Table 2), also known as AChE inhibitors (i.e., the test set).
The test set was compiled of organophosphates, carbamates and structurally non-classified pesticides. Regarding organophosphates, chlorpyrifos [37], DDVP [14], fenitrothion [38], azinphos-methyl [39], naled (dibrom) [40], TCVP [41], parathion [39], methyl parathion [42], diazinon [43], phosmet [14], azamethiphos [44], and terbufos [45], respectively, are used as insecticides. Moreover, glyphosate is a broad-spectrum systemic herbicide and crop desiccant [46], malathion is mostly used for mosquito eradication [47], parathion is used as an acaricide. Phosmet is mainly used on apple trees to control codling moth and on a numerous fruit crops, ornamentals, and wines for the control of aphids, suckers, mites, and fruit flies. Terbufos is a nematicide used on corn, sugar beets and grain sorghum, and it is extremely toxic to birds. On the other hand, carbamate pesticides, such as methomyl and methiocarb are used against insects; oxamyl is an extremely toxic pesticide to humans, fish, and birds [48,49,50], whereas methiocarb is additionally used as a bird repellent, acaricide and molluscicide. Among the structurally non-classified pesticides, 2,4-D is widely used as a herbicide, which mimics the action of the plant growth hormone auxin, resulting in the uncontrolled growth and eventual death in susceptible plants [14]; DDT is originally developed as an insecticide, used to control malaria and typhus, and as a contact poison against several arthropods [51]; DEET is the most common active ingredient in insect repellents for protection against mosquitoes, ticks, fleas, chiggers, leeches and many biting insects [52].
Among all of the listed pesticides (training and test set compounds), the pharmacology of organophosphorus (OP), compounds as hAChE inhibitors, is by far the most investigated; it involves the phosphorylation of catalytic Ser203 and the formation of stable inactive phosphyl-AChE adducts [36]. Moreover, methylcarbamate (MC) insecticides like carbofuran and similar compounds, occupy the catalytic triad area and act as reversible competitive inhibitors of ACh [73], elucidated by molecular docking (SB) studies on Nephotettix cincticeps AChE. Docking studies were also performed for imidacloprid analogues [74].
Therefore, to enrich the information related to the pharmacology of commonly used commercial pesticides, the targets of herein study, their mechanisms of action were either confirmed or nominated by means of reversible and reversible-like (for organophosphorus pesticides only, docking performed to obtain non-covalent conformation that is going to be transformed to the covalent one) cross-docking, as well as by single step molecular dynamics studies on Mus musculus and Homo sapiens AChE level. Alongside with the prediction of free energy of binding, generated bioactive conformations, were the starting points for the definition of the mechanism of action of 2-chloro-1,3,5-triazine-based pesticides (atrazine, simazine, and propazine, respectively), as the most toxic ones. DFT quantum mechanics mechanistic studies were supposed to shed further light on pesticide-AChE interactions for the purpose of designing pesticides with lower acute toxicity profiles.

2. Results and Discussion

2.1. The Conformational Analysis of Pesticides in Aqueous Solution

Conformational analysis (CA) studies are herein reported for simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] (Table 3), as pesticides with known crystal structures, as well as for atrazine and carbofuran, as pesticides that exert the highest acute toxicity (Table 4). To avoid redundancy, applied studies to the other pesticides are reported as Supplementary Material. The crystal structures of simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] (Table 3, Tables S1–S3) were used as starting conformational geometries for the CA. For pesticides without experimental data (Table 4, Tables S1–S3), starting conformations were built from scratch and energy minimized. For the sake of clarity, pesticides were divided into four groups based on their common structural patterns: 2-chloro-1,3,5-triazine derivatives: simazine, atrazine, and propazine (Table 3 and Table 4, Table S1, respectively); amide derivatives: monocrotophos, dimethoate, and carbofuran (Table 3 and Table 4, Table S2, respectively), carbaryl and tebufenozide (Table S2); 6-chloropyridine-3-yl)methanamine derivates: acetamiprid (Table 3, Table S3) and imidacloprid (Table S3); 1-(4-chlorophenyl)-3-methylurea condensates: diuron, monuron, and linuron (Table S3).
CA in explicit solvent may be used to predict the single pesticide stability prior to AChE binding. Herein, CAs were tentatively carried out by using five different force fields (FFs): MM3, AMBER94, MMFF, MMFFs, and OPLSAA [35]. Either for crystallized pesticides or pesticides with unknown experimental structures, the procedure was divided into the following steps: (1) conformational analysis; (2) determination of energy of the global minimum (Eglob_min), and (3) the LB superposition of structures to define the alignment accuracy. All steps were performed in order to obtain the best performing FF, either by means of capability to predict the free energy of global minimum or to superimpose molecules. The higher goal of CA was to answer the question: How to consider the FF of choice while treating pesticides: by predicted Eglob_min values or by alignment capability?
CA results for simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] are reported in Table 3. MM3 and AMBER94 could not be applied to all of the crystallized pesticides as the needed parameters for the optimization were not available. Possible workaround could be the utilization of Antechamber suite [75], by which the generic AMBER parameters could be calculated, but then a different protocol would have been developed instead of the one based on the well-known accredited Schrödinger software [76]. MMFF showed to be the only FF containing the parameterizations for all of the compounds; MMFFs failed in the optimization of simazine whereas OPLSAA was not applicable to monocrotophos and acetamiprid.
The absolute values of Eglob_min are not necessarily good measures of FF quality [77]. FFs and are usually parameterized to give good values of relative energies among different conformations of the same molecule [77], so a lower value of energy for a global minimum does not necessarily mean that a given FF would perform better [77]. The energy differences between the conformations are simple reliable indicators of conformational analysis algorithm ability than the absolute energy values [77].
Therefore, the evaluation of conformation reproducibility was achieved by means of calculated minima root mean square deviation (RMSD) values [78] upon conformations’ superimposition (Table 3 and Table 4, Tables S1–S3). Comparison of the experimentally available crystal structures conformations with those obtained by CA revealed, a high level of performance by MMFF (displaying RMSD values lower than 1 Å for all of the four considered pesticides) [78]. That was somehow expected, given that MMFF-based optimization was applicable to all of the listed pesticides. Regarding the MMFFs and OPLSAA, these FFs only succeeded in reproducing the crystal structure of acetamiprid.

2.2. The Prediction of Pesticides Structures in Aqueous Solution

The above assessed CA procedure was further applied to the training set pesticides with an unknown crystal structures. The results of optimization for the pesticides with the highest acute toxicity (i.e., the lowest LD50 values), atrazine and carbofuran, using the best performing FFs, are reported in Table 4, while the comparison of all utilized FFs for atrazine, carbofuran, and the remaining compounds are presented in Tables S1–S3. As expected, AMBER94 FF, as implemented in Schrödinger suite [76] was not applicable at all, MM3, MMFF and MMFFs FFs gave similar global minima for the majority of the listed pesticides, while OPLSAA was not applicable for all compounds. A significant conformational alignment accuracy (RMSD < 1 Å) was achieved for atrazine (MM3/MMFF, MM3/MMFFs, and MMFF/MMFFs combinations, Table 2, Table S1) and carbofuran, (MMFF/MMFFs, MMFF/OPLSAA, and MMFFs/OPLSAA combinations, Table 3, Table S2). Summarizing the CA results, conclusion can be made, despite the limited number of training set compounds, MMFF is universally applicable FF for either Eglob_min calculation or RMSD-based alignment accuracy assessment and can be used in any further study, describing the pesticides conformational analysis.

2.3. The Prediction of Externally Evaluated Pesticides Structures in Aqueous Solution

MMFF was thereafter used to predict the global minimum conformation of test set pesticides (Table 2) as well. The results of the optimisation are presented as Supplementary Materials (Table S4). Global minimum conformations were generated for by the means of the identical CA setup. Generated conformations were then used for the purposes of training set results external validation.

2.4. The Pesticides Acute Toxicity Anticipation through the Pre-Bound Conformations

The lowest energy conformers of the 2-chloro-1,3,5-triazine group were of high stability in the aquous solution, according to the values obtained by MMFF (Table 3 and Table 4, Table S1). Being chemically stable, pesticides may in perspective saturate the AChE active site and low values of global minimum (Eglob_min) represent the indicators of high acute toxicity (i.e., the low LD50 values) against Mus musculus. For either atrazine, propazine, or simazine, low Eglob_min values are indeed the precondition of high acute toxicity (Table 1: LD50 values equal to 0.85, 3.18, and 5 mg/kg of body weight, respectively, Table 3 and Table 4, Table S1: Eglob_min values equal to −1007.32, −1022.18, and −992.47 kJ/mol, respectively). For much of the lasting pesticides (monocrotophos and dimethoate: Table 1 vs. Table 3 and Table S2; carbaryl and tebufenozide: Table 1 vs. Table S2; imidacloprid: Table 1 vs. Table 3 and Table S3, and acetamiprid: Table 1 vs. Table S3), it appeared that high to medium acute toxicities mostly correlate with low to medium Eglob_min values.
The above paradigm was confirmed for carbaryl and tebufenozide, where the high Eglob_min values are a precondition for low acute toxicity (Table 1: LD50 values equal to 100 and 5000 mg/kg of body weight, respectively, Table S2: Eglob_min values equal to 20.43 and 467.38 kJ/mol, respectively), as well as for monocrotophos and dimethoate, whose medium LD50 values (Table 1: LD50 equal to 14 and 60 mg/kg of body weight, respectively) are associated with medium Eglob_min values (Table 3: Eglob_min values equal to −307.36 and −383.68 kJ/mol, respectively). On the other hand, the high acute toxicity of carbofuran against Mus musculus (Table 1: LD50 = 2 mg/kg of body weight) disagreed with calculated high Eglob_min values (Table 4: Eglob_min value equal to 20.30 kJ/mol). Interestingly, a thermodynamic parameter calculated for carbofuran suggests some instability in aqueous solution; according to the previous observations, a high dosage of pesticide would be required to induce mortality in Mus musculus. Therefore, other parameters need to be considered, to explain why carbofuran acts as an outlier.
Among the last group of pesticides, acetamiprid (Table 1 vs. Table 3 and Table S3) and imidacloprid (Table 1 vs. Table S3) also follow the archetype that relatively low Eglob_min values are a precondition for medium LD50 dosages. Linuron, however, behaves oppositely to 2-chloro-1,3,5-triazines group inasmuch as with its low Eglob_min values there is a high LD50 value (Table 1: LD50 = 2400 mg/kg of body weight, Table S3: Eglob_min value equal to −992.47 kJ/mol). Linuron’s structural analogues diuron and monuron, despite the fact they are slightly less stable (higher Eglob_min values, Table S3: Eglob_min values equal to −142.96 and −197.63 kJ/mol, respectively) need to be administered at high dosages to induce mice mortality (Table 1, LD50 = 500 and 1700 mg/kg of body weight, respectively).
In order to enumerate the Mus musculus LD50/Eglob_min interconnection, simple linear regression (as implemented in LibreOffice Calc) was analysed considering pLD50 (calculated by converting LD50 from mg/kg of body weight to −log(mol/kg of body weight)), as dependent variable, and Eglob_min values as calculated by MMFF (Table S5), as descriptive independent variables. The use of all the data (model 0) did not lead to a significant statistical regression. However, the exclusion of carbofuran and linuron as highlighted outliers (those with the worse pLD50 fitted values), led to a monoparametric QSAR model 1 (Equation (2)):
  p L D 50   = 0.0021 E glob _ min 2.82   ( r 2 = 0.78 ) ,
Additional exclusion of diuron and monuron further increased the statistical robustness leading to a QSAR model 2 (Tables S5 and S6, Figure 2) with more than 90% of explained variance (90.5%) (Equation (3)):
  p L D 50   = 0.0019 E glob _ min 3.05   ( r 2 = 0.905 ) ,
The obtained model 2 estimated the pLD50 values of carbofuran, linuron, diuron, and monuron, as in situ generated test set, with good accuracy (Figure 2a). On the other hand, model 2 was also used to predict the Mus musculus pLD50 values of the test set pesticides (Table 2 and Table S6, Figure 2b). However, despite the good statistical potential of model 2 (i.e., the explained variance value), the Eglob_min values cannot be considered as valuable indicators to anticipate the Mus musculus pLD50 values and for the acute toxicity prediction.
Moreover, models 1 and 2 cannot be universally applied for the prediction of the Homo sapiens pLD50 as the Eglob_min values are common to both model organisms. Therefore, different models are required for Homo sapiens as model organism. In that sense, additional models were developed from calculated Homo sapiens pLD50 values (Equation (1)). As expected, the incorporation of all the training set, pesticides did not lead to any statistically significant model. As above, exclusion of carbofuran and linuron led to the model 3 (Equation (4)):
            p L D 50   = 0.0021 E glob _ min + 3.01   ( r 2 = 0.783 ) ,
and further exclusion of diuron and monuron led to model 4 (Equation (5), Tables S7 and S8, Figure 2c,d):
    p L D 50   = 0.0019 E glob _ min + 4.15   ( r 2 = 0.91 ) ,
Like for Mus musculus and Homo sapiens, the pLD50 values were not fully satisfactory correlated to the global minima energies, leading to the conclusion that Eglob_min values are not adequate indicators to anticipate the pLD50 against any model organism, the level was changed from LB to SB, directing that the free energies of binding could be considered as appropriate acute toxicity indicators. Free energies of binding can be calculated separately for each pesticide and model organism, by performing molecular docking studies, either for Mus musculus AChE or Homo sapiens AChE.

2.5. The Structure-Based Studies

To get more insights into the mechanistic profile of understudied pesticides, their interactions with either mAChE and hAChE, as molecular targets, were also considered. Sequence similarity value between mAChE and hAChE was found to be 88.44% (Figure S1), rising to the 100% identity between [the active sites, indicating that similar pharmacodynamics profile should be expected with either mAChE or hAChE by the same molecule. Therefore, the application of SB tools could represent the computational approach to predict the biochemical mechanism and hence the level of acute toxicity. Two SB methodologies were chosen: molecular docking and molecular dynamics. Docking and dynamics studies were initially performed to anticipate the bioactive conformations of pesticides with known acute toxicity and unknown binding mode into the AChE, after which the obtained free energies of binding were correlated to the acute toxicities. Despite the fact there are no co-crystalized AChEs in complex with pesticides as either non-covalent or covalent inhibitors, it was arbitrarily decided to simulate reversible docking, for pure reversible AChE inhibitors, and reversible-like docking, for covalent AChE inhibitors, followed by the appropriate, one step molecular dynamics simulations. Reversible-like docking, was supposed to obtain the pre-covalent conformation that will lead to the AChE covalent modification by literature-claimed [73] covalent inhibitors. The ultimate task of SB studies was to reveal the free energies of binding and to correlate them to the pesticides acute toxicities.

2.5.1. The Alignment Rules Assessment

In order to learn how to perform the SB alignment of targeted pesticides within the AChE active sites of either mouse or human molecular target, the 3-D coordinates of reversible AChE inhibitors were taken from their corresponding minimized complexes and used to gain the knowledge how to reproduce their co-crystallized conformations by means of the SB alignment assessment, as well as to determine the best SB molecular docking methodology. Within the SB alignment assessment, AutoDock [79], AutoDock Vina [80] (in the future text just Vina) and DOCK [81] were evaluated to select the best docking algorithm/scoring function and to perform the SB positioning of pesticides. The SB alignment assessment procedure was performed with the standard four step protocol [78]: Experimental Conformation Re-Docking (ECRD), Randomized Conformation Re-Docking (RCRD), Experimental Conformation Cross-Docking (ECCD), and Randomized Conformation Cross-Docking (RCCD).
Thus, within the ECRD stage, Vina and DOCK possessed the highest ability to reproduce the experimental conformations of Mus musculus and Homo sapiens AChE inhibitors. Even 75% of the available inhibitors were correctly reproduced by Vina, whereas DOCK was slightly less potent, reproducing the 50% of structures in either rigid or flexible manner, respectively (Tables S9 and S10). Similarly, Vina was of the highest potency in the initial reproduction of Mus musculus AChE inhibitors with the docking accuracy (DA) of 55.55% (Table S9). DOCK, however, completely failed in the reproduction the experimental poses of Mus musculus AChE inhibitors (Table S9). In the initial process, as well as in all other difficulty levels, AutoDock algorithm failed in the experimental conformations reposition.
Similar potency for Vina and DOCK was observed within the RCRD stage in the process of human inhibitors reproduction (Table S10, Vina DA equal to 68.75%, DOCK DA equal to 56.25%, respectively). While processing mouse inhibitors (Table S9), Vina re-positioned the modelled inhibitors structures with the high efficiency (DAs equal to 57.78%), while DOCK’s potency dropped below 30%.
The level of difficulty was increased during the ECCD stage where either Vina or DOCK retained the high level of accuracy while cross-aligning of human inhibitors (Table S10). As for the other model organism, either Vina or DOCK were inaccurate while matching the co-crystallized inhibitors with diverse mouse enzymes (Table S9).
The decisive difference in favour of Vina, by means of the general docking accuracy, was noticed within the RCCD stage, as the most difficult one. The cross-alignment of randomized conformations is considered as the most important validation stage for any docking program as it represents the ability of a program to position a random ligand conformation in a non-native environment, i.e., into its molecular target whose active site amino acids suffered induced fit conformational change due to the presence of structurally different inhibitors. Within the RCCD stage, Vina was far superior in comparison with other programs, by means of The SB alignment accuracy, and, therefore, can be considered as a valuable tool for the prediction of bioactive conformations of acetylcholine esterase inhibitors with the unknown binding mode.

2.5.2. The Cross-Docking and Molecular Dynamics Studies of ACh

Both the training set or the test set pesticides were submitted to the herein Vina-based proposed SB alignment protocols, that have been validated by reproducing the bioactive conformations of PDB-available inhibitors in complex with AChE (Supplementary Materials Tables S9 and S10, and Figure S2) and applied in situ in predicting the bioactive conformations of ACh (Figure 3).
As no experimental structures of ACh in complex with either wild-type mAChE or hAChE are available in the Protein Data Bank, its bioactive conformation within both enzymes was predicted by means of Vina, as the best performing SB alignment assessment tool (Figure 3a,b, see Supplementary Materials Alignment Rules Assessment section, Tables S9 and S10). The stability of either ACh-mAChE or ACh-hAChE complex was confirmed using one step molecular dynamics (MD) runs, upon the inspection of Cα RMSD values (1.2 ns simulation length, Figures S7a and S8a), RMSF values (Figures S7b and S8b), and trajectories energy terms (data not shown but available from the authors upon request). The short-termed MD protocol was sufficient to produce stable complexes. The two ACh poses (in mAChE or hAChE) were conformationally very close, displaying a RMSD value of 0.043 Å.
As expected, the positive quaternary ACh amine group was located in the anionic sub-site, where the orientation of particular group is supported by strong hydrophobic interactions with a m(h)Trp86 methyl group, as well as with m(h)Tyr124, m(h)Tyr337, and m(h)Tyr341 aromatic moieties (Figures S7c and S8c). The positively charged amine was electrostatically attracted by the electronic cloud of the m(h)Trp86 aromatic group. The acetyl group occupied the esterase sub-site: the carbonyl oxygen was involved in the electrostatic interactions with m(h)Tyr337 side chain hydroxyl group and m(h)His447 imidazole ring, while the carbonyl oxygen was oriented towards the side chain of m(h)Tyr133. Regarding the mouse enzyme, during the MD run, the carbonyl moiety was observed to be involved in the hydrogen bonding with mTyr133 for 29% of simulation time (Figure S8c, dHB = 1.943–2.337 Å) whilst for human enzyme the carbonyl group was mainly involved in the electrostatic interactions with hSer203 and hHis447 (Figure S8c). The analysis of docking and MD results agreed with the hydrolysis reaction pathway of ACh (Figure 1a) that was also externally confirmed by QM/MM studies elsewhere [82]. In fact, the ACh molecule had been approaching m(h)Ser203 during the simulation for the 57% of the time, where the range of distances between the acetyl group and the catalytic residue was 1.473–1.894 Å (Figures S7c and S8c).

2.5.3. The Cross-Docking and Molecular Dynamics Studies of Pesticides 2-chloro-1,3,5-triazine Group

Atrazine’s, simazine’s and propazine’s acute toxicity (Table 1) can be explained by means of their docked conformations, within either mAChE or hAChE (Figure 4 and Figure 5, respectively). The graphical interpretation of the results of the SB and MD studies are herein shown for atrazine (Figure 4a and Figure 5a, Figure S9 and Figure S10, respectively) and simazine (Figure 4b and Figure 5b, Figure S13 and Figure S14, respectively), whereas the corresponding illustrations considering the molecular modelling studies of propazine are reported as Supplementary Materials (Figure S3a, Figure S5a, Figure S11, Figure S12, respectively).
The best-docked poses of the triazine-based pesticides were found to occupy the ACh binding cavities of both the mAChE and hAChE. The main differences were not related to the predicted pesticides bioactive conformations, but were mainly ascribed to enzymes active site residue Thr83: the side chain methyl group of mAChE Thr83 pointed toward the anionic sub-site, whereas in the human enzyme the corresponding group is flipped out vertically from the active site. The latter caused the differences in the stabilization period during the MD runs: atrazine-mAChE, propazine-mAChE, and simazine-mAChE complexes were respectively of high stability after 0.100 ns, 0.120 and 0.610 ns (Figures S9a, S11a and S13a), whereas the complexes of human analogues equilibrated approximately 20 ps afterwards (Figures S10a, S12a and S14a).
Nevertheless, within either mice or human bioactive conformations, the 2-chloro-1,3,5-triazines were parallel to m(h)Trp86 and m(h)Tyr337, trapped in the hydrophobic cage between two residues. The individual contribution of these interactions was respectively −37.34 and −19.9 kcal/mol, according to the Free-Energy Decomposition Analysis, FEDA (Table S11). Atrazine bioactive conformations were nearly identical in both enzymes. On the other hand, comparing simazine docked structures (Figure 4b and Figure 5b), the one found in hAChE (Figure 5b) was in-plane rotated for about 30°. The triazine nitrogen atom, p-positioned in relation to chlorine, was electrostatically stabilized by both the hydroxyl group of m(h)Tyr337 and the indole nitrogen of m(h)Trp86. The individual contribution of these interactions was respectively −34.26 and −15.56 kcal/mol, according to FEDA (Table S11). Within the mAChE, the m-N’ atoms, holding the isopropyl function of atrazine (Figure 4a), and one of the ethyl groups within simazine (Figure 4b), were oriented toward the mThr83 methyl group and mTyr341 phenolic moiety. All the described interactions were confirmed during the MD simulation (Figure S9c and Figure S13c, respectively). The atrazine m-N’ atom was involved in a strong hydrogen bonding (HB) (dHB = 2.253–2.976 Å) with the hThr83 hydroxyl group (Figure S10c) in 56% of MD run time. A similar interaction was not detected during the atrazine-mAChE complex MD run (Figure S9c), suggesting that the hThr83 conformational flip and the existence of hThr86-m-N′ HB may even increase the acute toxicity of atrazine to humans. Regarding the simazine, as a consequence of the in-plane rotation, its m-N’ atom created HB with hTyr341 (dHB = 2.341–3.127 Å, 76% of the time of simulation, Figure S14c). Furthermore, the functional groups that substitute the remaining m-N″ atoms were directed towards the hTyr133, hSer203, and hHis447, and HB-attracted by hHis447 (simazine established the strongest HB, dHB = 2.379–3.242 Å, 94% of MD time, Figure S14c). Since the precise hydrogen bonds were not observed during the MD simulation within the mAChE, they may represent some indicators that the atrazine and simazine could exert even higher acute toxicity against humans than against mice. Moreover, the described docking/MD interactions were all confirmed by the FEDA application (Table S11).
Regarding propazine, the m-N′ atoms, holding the isopropyl units or propazin (Figure S5a) were oriented within the mAChE towards the side chain methyl group of mThr83 and a phenolic moiety of mTyr341. Similarly to atrazine and simazine, propazine may be even more toxic to humans than to mice, as it exerts the toxicity against humans by establishing the hydrogen bonds with hTyr341 (dHB = 2.286–3.114 Å) and hHis447 (dHB = 2.121–3.331 Å) in the 45% and the 95% of time, respectively (Figure S12c). Attracted to mThr83, propazine was limited with the ability to make hydrogen bonds with mTyr341 or mHis447, as confirmed by MD studies (Figure S11c).

2.5.4. The Cross-Docking and Molecular Dynamics Studies of Pesticides Amide Group

Among the amide-based pesticides, carbofuran and monocrotophos were of the highest acute toxicity (Table 1). Therefore, the results of carbofuran (Figure 4c and Figure 5c, Figure S15 and Figure S16, respectively) and monocrotophos (Figure 4d and Figure 5d, Figures S17 and Figure S18, respectively) docking and MD studies are herein reported and discussed, whereas the corresponding graphical interpretations for lower toxic dimethoate (Figure 4e, Figure 5e, Figure S19, Figure S20, respectively), carbaryl and tebufenozide are reported as a Supplementary Materials (Figure S4, Figure S5, Figures S21–S24, respectively).
Carbofuran was the most toxic against Mus musculus (Table 1). According to the literature, this particular pesticide is, as a carbamate derivative, proclaimed to be the reversible AChE inhibitor [73]. The results of herein conducted SB studies are in full agreement with the proposed/confirmed reversible inhibition of AChE by carbofuran. Hence, within either mAChE or hAChE active site, the compound’s methyl methylcarbamate scaffold contributed as follows: the N-methyl group was positioned between the methyl group of m(h)Thr83 and m(h)Trp86 (Figure 4c and Figure 5c), a carbonyl group interfered with m(h)Tyr341 while the ether link was electrostatically stabilized by hydroxyl groups of m(h)Tyr124, m(h)Tyr337, and m(h)Tyr341, respectively. However, during the MD simulations, methylcarbamate displayed rotation towards the m(h)Tyr124 (Figures S15c and S16c) upon which the amide carbonyl group was involved in strong HB with m(h)Tyr124 for the 40% of time (dHB = 2.478−2.783 Å, Figure S15c and Figure S16c, respectively), indicating a possible reason for carbofuran acute toxicity. The other HB was formed between the secondary amine and m(h)Trp86 (dHB = 2.784−3.226 Å, Figure S15c and Figure S16c, respectively). Furthermore, the 2,2-dimethyl-2,3-dihydrobenzofuran part of carbofuran was located beneath the m(h)Trp86 and established hydrophobic interactions with m(h)Tyr337 due to the π-π stacking, as well as the electrostatic interactions with m(h)Tyr124 via 2,2-dimethyltetrahydrofuran oxygen atom. Moreover, all the described docking/MD interactions were confirmed upon conducting the FEDA procedure (Table S11).
The covalent mode of inhibition [73] and the high acute toxicity of organophosphorus-based pesticides monocrotophos and dimethoate against Mus musculus and Homo sapiens were primarily conditioned by the interactions of their amide groups with AChE. Hence, in the best docking pose of monocrotophos within mAChE, the amide carbonyl was oriented toward the hydroxyl group of mTyr341 (Figure 4d). Such interaction was a result of hydrophobic stabilization between the pesticide’s allyl methyl group and the mTyr341 phenyl ring, as well of the N-methyl group and mThr83. Still, within the human enzyme (Figure 5d), corresponding amide carbonyl group approached the hThr83 hydroxyl group, confirming that the orientation of hThr83 makes a difference in SB alignment in hAChE. The consequent important monocrotophos interactions were those generated by the dimethylated phosphoric acid. Thus, the phosphoric acid P=O group established electrostatic interferences with mAChE Tyr124 (Figure 4d); the precise oxygen atom served as an HB-acceptor to mTyr124 during the MD simulation (dHB = 2.227–2.946 Å, Figure S17c). For the purpose of hAChE inhibition, the corresponding monocrotophos carbonyl group was, due to the strong hydrophobic attraction between the phosphoric acid methoxy groups and hTrp86, only partially directed towards the hTyr124 (Figure 5d). However, monocrotophos adapted its conformation/binding upon 0.718 ns of MD simulations with hAChE (Figure S18a), and, contrarily to mAChE, hTyr124 was involved in 96% of the time in HB with the amide carbonyl (dHB = 1.984–2.117 Å, Figure S18c). Moreover, one of the phosphoric acid oxygen atoms made HB with hTyr337 (35% of simulation time, dHB = 2.449–2.688 Å, Figure S18c). Still, inhibiting mAChE or hAChE, the monocrotophos phosphoric acid moiety remained situated in the proximity of m(h)Ser203 and m(h)His447 (the average distances: 1.25–1.783 Å for mAChE, 1.375–1.651 Å for hAChE, Figure S17c and Figure S18c, respectively), suggesting that the covalent modification of mSer203 or hSer203 will occur in vivo [74]. Moreover, all the described docking/MD interactions were confirmed upon conducting the FEDA (Table S11).
The monocrotophos structural analogue, dimethoate, is, as Mus musculus and Homo sapiens AChE inhibitor, limited with the capacity to interfere with m(h)Tyr341 by means of steric interactions, inasmuch as it does not contain the allyl methyl group (Figure 4e and Figure 5e, respectively). Hence, although dimethoate is in both active sites superimposed to monocrotophos, by means of backbone orientation (Figure 4d and Figure 5d, respectively), its amide carbonyl is oriented towards m(h)Tyr124 and m(h)Tyr337. The bioisosteric replacement of the phosphoric acid carbonyl group with the thione one influenced that the sulfur atom is within the mouse enzyme directed to the hydroxyl group of m(h)Tyr337, while within the human one the thione moiety is in reach of the m(h)Ser203 and m(h)His447.
However, surprisingly, there was a little consistency between the best-docked structure of dimethoate within the mouse active site (Figure 4e) and the molecular dynamics solutions (Figure S19). Thus, pesticide’s amide carbonyl established no interactions, while the amide nitrogen was periodically involved in hydrogen bonding with mTyr341 time (dHB = 2.971–3.557 Å, Figure S19c). Thione functional group was HB-acceptor to mHis447 (25% of the time, dHB = 2.894–3.431 Å), while one of the methoxy functions was in the HB-vicinity to mGly121 (27% of the time, dHB = 2.723–3.338 Å). Those interactions allegedly provide the opportunity for mSer203 to perform the electrophilic attack on pesticide (Figure 6f). On the other side, there was a high level of agreement between the dimethoate rigid docking (Figure 5e) and molecular dynamics within the human enzyme, by means of the amide carbonyl group interactions, given that this pharmacophore was for the 95% of time included in hydrogen bonding with hTyr124 (dHB = 1.643–2.128 Å, Figure S20c). Moreover, one of the phosphoric acid oxygens was, as for monocrotophos, used to make a hydrogen bond with hTyr337 (35% of simulation time, dHB = 2.449–2.688 Å, Figure S20c). In general, considering the overall orientation of the molecule, dimethoate would most likely behave as a covalent inhibitor of human AChE (Figure 6f), as well.
Carbaryl reached the comparable bioactive conformation to carbofuran (Figures S3b and S5b), within both species AChE active sites, with the difference that the 2,2-dimethyltetrahydrofuran scaffold is in the structure of particular pesticide replaced with another benzene ring. Such replacement has a direct consequence in a diminished toxicity of carbaryl in comparison with carbofuran, against either mice or humans, since carbaryl cannot establish additional electrostatic interactions with m(h)Tyr124.
Regarding carbaryl, it is proclaimed (Table 1) that this pesticide is more than fifty times less toxic than carbofuran when administered orally to Mus musculus, and MD revealed the reason. Hence, within the mAChE, carbaryl’s amide carbonyl group forms a hydrogen bond, not with the mTyr124, but with the HB-network made of ordered water molecules and mAsn87 (Figure S21c). The precise function also establishes the pure HB with mSer125 (75% of simulation time, dHB = 2.335–2.412 Å, Figure S21c). However, according to the MD study performed on the carbaryl-hAChE complex, there is a 36% of probability that particular carbonyl group is oriented towards HB with hTyr124 (dHB = 2.374–2.894 Å, Figure S22c), thus anticipating the high toxicity of this pesticide to humans.
The lowest active amide-based pesticide is tebufenozide (Table 1). Almost symmetrical as a compound, this pesticide can be divided in two scaffolds by which it establishes the interactions within the active site of either mAChE or hAChE: the 4-ethylbenzamide and the N-tert-butyl-3,5-dimethylbenzamide (Figure S3c and Figure S5c, respectively). The RMDS value between the mouse and human docked conformations amounts 0.472 Å, implying the high similarity between the obtained conformations. Thus, the ethylbenzene scaffold of 4-ethylbenzamide was oriented parallel to benzene core of m(h)Tyr337, establishing hydrophobic interactions. The carbonyl group that completes the 4-ethylbenzamide was electrostatically attracted by three active site residues: m(h)Tyr124, m(h)Tyr337 and m(h)Tyr341. The hydrophobic part of m(h)Tyr124 additionally stabilized the tert-butyl group, incorporated in the second part of pesticide. The carbonyl group belonging to the tert-butyl-3,5-dimethylbenzamide moiety was geared towards the amide group of m(h)Gln71 and the hydroxyl group of m(h)Tyr124, for the purpose of forming the electrostatic interactions. The remaining part of the distinct subsection, in the form of 3,5-dimethylbenzene, is stabilized with the hydrophobic parts of m(h)Trp86 and m(h)Phe295. Similarly to other amide-based pesticides, this pesticide exerts the toxicity after creating the strong hydrogen bond in-between the tert-butyl-3,5-dimethylbenzamide moiety and either mTyr124 or hTyr124 (95% of simulation time, dHB = 1.273–1.588 Å, Figure S23c; 94% of simulation time, dHB = 1.268–1.575 Å, Figure S24c). Beside those beneficial interferences, the pesticide owns its stability and consequent toxicity to numerous hydrophobic interactions (Figure S23c and Figure S24c, respectively).

2.5.5. The Cross-Docking and Molecular Dynamics Studies of Pesticides (6-Chloropyridin-3-Yl)Methanamine and 1-(4-Chlorophenyl)-3-methylurea Groups

The pesticides comprising the two remaining groups exerted the lowest acute toxicities to Mus musculus (Table 1). The results of SB and MD studies are herein shown for imidacloprid (Figure 4f and Figure 5f, Figure S25 and Figure S26, respectively), while the corresponding exemplifications considering the molecular modelling studies for lower toxic acetamiprid, diuron, monuron, and linuron are available as Supplementary Material (Figures S3–S6 and Figures S27–S34, respectively).
Regarding the best docking poses of imidacloprid within either mAChE or hAChE, these indicated that the 2-chloro-5-methylpyridine moiety was sterically stabilized by m(h)Trp86 and m(h)Tyr337 (Figure 4f and Figure 5f, respectively). The difference in the alignment occurred with the N-(imidazolidine-2-ylidene)nitramide scaffold, whose nitro group was oriented towards mSer203, while the same function in human enzyme was in the proximity of hTyr124. Even though the imidacloprid-mAChE complex was more stable in comparison to the human one (Figure S25a), interactions generated within the complex were inconsistent to those that imidacloprid established with hAChE. Thus, while interacting with mAChE (Figure S25c), the nitro group was involved in the hydrogen bonding with mGly122 (64% of MD time, dHB = 2.483–2.788 Å) and mAla204 (59% of MD time, dHB = 2.771–3.132 Å). The corresponding functional group in hAChE (Figure S26c) was a HB acceptor for hTyr124 (82% of MD time, dHB = 1.664–1.738 Å) and hTyr337 (68% of MD time, dHB = 1.825–1.931 Å), and further confirmed those residues as essential for pesticides acute toxicity in humans. Moreover, all the described docking/MD interactions were confirmed by FEDA (Table S11).
The important interactions anticipated by the molecular docking and molecular dynamics (Figure 6), provided the basis for the understanding of the general mechanism of acute toxicity on the reversible and reversible-like inhibitors level. Even though the covalent AChE inhibitors are highly efficient as pesticides, the future pesticides production must be based on the reversible-like interactions, avoiding the pesticide-AChE complex regeneration or ageing [73]. The obtained reversible-like docking and MD poses for covalent AChE inhibitors clearly confirmed the potency of selected computational tools (Vina/Desmond programmes pair) to predict the covalent binding mode of any compound in future considered as covalent AChE inhibitor (i.e., the covalent pesticide). Therefore, there was no point to conduct the covalent docking that would be reduced just to a conformational search of the bonded molecules. Hereby applied docking/MD application was also aimed to assess the docking/MD software real utility to investigate the covalent AChE inhibitors. It must be emphasized that any inhibitor, even covalent, must undergo equilibrium between the free and the bonded stage. Therefore, the reversible-like (pre-covalent) pose is of extreme utility to predict the ligand conformation that will proceed to make a covalent bond. This information could be used in future works to design new reversible inhibitors (pesticide).
Considering the docking poses of acetamiprid (Figure S3d and Figure S4d, respectively), pesticide’s 2-chloro-5-methylpyridine moieties were aligned to each other in both enzymes and located in-between the m(h)Trp86 and m(h)Tyr337. The toxicity of particular pesticide may be additionally seen through the interactions of (E)-N-ethyl-N′-ethynyl-N-methylacetimidamide function, linked directly to the 2-chloropyridine ring and located in the spatial pocket made of m(h)Trp86, m(h)Gly121, m(h)Tyr133, m(h)Ser203, and m(h)His 447. However, according to the molecular dynamics results (Figure S27c and Figure S28c, respectively), this part of the pesticide has not achieved significant interactions with the specified site of amino acids, and the low toxicity of acetamiprid may be attributed only to the hydrophobic interactions established between the 2-chloro-5-methylpyridine and m(h)Trp86 and m(h)Tyr337.
Concerning the remaining compounds, diuron, linuron, and monuron, the basic scaffold they share is the 1-(4-chlorophenyl)-3-methylurea. In the structure of diuron, the terminal amine is dimethyl substituted while the aromatic moiety bears additional m-Cl substituent. Matching the conformations of diuron generated in Mus musculus and Homo sapiens AChE (Figure S4a and Figure S6a, respectively), conclusion can be made that the structures were almost perpendicular to each other, by means of aromatic rings positioning: the aromatic scaffold in mouse enzyme active site was perpendicular to mTrp86 and mTyr337, while in the human environment this scaffold adopted parallel orientation related to hTrp86 and hTyr337 (analogously to previous compounds). Considering the mAChE, the only reasonable explanation for such alignment of diuron is the interaction of mThr83 methyl group with the pesticide aromatic scaffold; the mThr83 side chain was positioned inside the enzyme active site and not out of it, like in humans (Figure S4a). Because of the different alignment, the carbonyl group of diuron′s 1,1-dimethylurea scaffold was oriented towards the mTyr124, mSer203, mTyr337, and mHis447 (Figure S4a). In the human environment (Figure S6a), the exact carbonyl group was directed towards the hTrp86 and hTyr133. The MD simulation of mouse complex confirmed the weak hydrogen bonding between the urea carbonyl group and mTyr124 (31% of simulation time, dHB = 2.872–3.132 Å, Figure S29c), while the corresponding interaction could not be established in human surroundings (Figure S30c). The low intensity of precise hydrogen bond may be the reason of low toxicity of diuron, latter on concluded for linuron, and monuron, too.
Speaking of the remaining derivatives, the aromatic scaffolds of monuron (Figure S4b and Figure S6b, respectively), and linuron (Figure S4c and Figure S6c, respectively) were in a sandwich between the m(h)Trp86 and m(h)Tyr337, in either mouse or human enzyme active site. The monuron carbonyl group was oriented towards the m(h)Tyr133 and m(h)Ser203, whereas the particular group in linuron is headed to m(h)Tyr 124, m(h)Tyr133 and m(h)Ser203. The nature of monuron and linuron carbonyl groups interactions with specified amino acid had been revealed after the MD studies. Thus, for either mouse or human AChE-based bioactive conformations, the carbonyl group of monuron interfered with m(h)Tyr124 in a manner of moderate strength hydrogen bonding (within the mouse enzyme active site during the 62% of simulation time, dHB = 2.548–3.012 Å, Figure S31c; within the human enzyme active site during the 56% of simulation time via the ordered water molecule, dHB = 2.936–3.326 Å, Figure S32c). Similar interactions were observed for linuron as well (within the mouse enzyme active site during the 56% of simulation time via the ordered water molecule, dHB = 2.985–3.455 Å, Figure S33c; within the human enzyme active site during the 90% of simulation time, dHB = 2.847–3.278 Å, Figure S34c).

2.5.6. The Cross-Docking of Externally Evaluated Pesticides

The structure-based studies by means of reversible and reversible-like bioactive conformations generation, recommended Vina as a tool for pesticides cheminformatics treating. Therefore, Vina was used to predict the global minimum conformation of test set pesticides (Table 2) as well.
The results of molecular docking are presented as a Supplementary Materials (Figures S35–S42). The generated conformations were used further on for the external validation of the training set results.

2.6. The Pesticides Acute Toxicity Anticipation through the Bioactive Conformations

Upon conducting the reversible and reversible-like molecular docking and molecular dynamics studies, for either training or test set pesticides, the acute toxicity was considered with regard to the established bioactive conformations. Therefore, similarly as above reported for global minima energies, monoparametric QSAR model was derived to initially describe the acute toxicity against Mus musculus by enumerating the acute toxicities through the free energies of formation of reversible and reversible-like binding for all the training set pesticides (model 5, Equation (6), Table S12):
        p L D 50   = 0.7885 Δ G binding 2.44   ( r 2 = 0.95 ) ,
A satisfactory model was obtained by means of the whole training set and was used to recalculate the pLD50 values for acute toxicity against Mus musculus, for both the training and test sets (Figure 7a,b, Tables S12 and S13). Afterwards, the model was used to predict the pLD50 values for acute toxicity against humans on the basis of ∆Gbinding values for Homo sapiens AChE (Figure 7c,d, Tables S14 and S15).
Thus, the QSAR model 5 described the acute toxicities against Mus musculus, by means of the free energies of binding of reversible and reversible-like bioactive conformations for either non-covalent or covalent pesticides, respectively. The model is also proven to be a tool for the anticipation of acute toxicities against humans. In the end, the conduced structure-based studies unequivocally supported pLD50/∆Gbinding interrelation obtained by QSAR model 5. In that sense, the SB approach, conducted through Vina as the cheminformatics engine, was proven to be the correct one for the anticipation of acute toxicities of targeted pesticides. Consequently, the Vina-based ∆Gbinding values can be considered in future design as valuable indicators to anticipate the pesticides LD50 values against any model organism.

2.7. The Quantum Mechanical Studies of Acetylcholinesterase Inhibition by 2-Chloro-1,3,5-triazine-based Pesticides

The performed molecular docking and molecular dynamics studies agreed with the previously reported mode of action of amide-based, 6-chloropyridine-3-yl)methanamine-based, and 1-(4-chlorophenyl)-3-methylurea-based pesticides [73]. The quantum mechanical (QM) calculations were only performed on hAChE, with the aim to reveal the pharmacology of 2-chloro-1,3,5-triazine-based pesticides and the origin for the acute toxicity, given that their pharmacology is the least understood.
As molecular docking and molecular dynamics studies indicated the formation of hydrogen bonds network between the atrazine, propazine, and simazine and hAChE, the HB-bearing pesticides-hAChE MD poses were extracted from pesticides-hAChE complexes and subjected to quantum mechanical (QM) DFT mechanistic studies. Thus, the DFT-based calculations were performed on two different levels. The first level included the extraction of MD geometries where favourable hydrogen bonds between the 2-chloro-1,3,5-triazine-based pesticides and hThr83, hTyr341, hTyr124, and hHis447, respectively, were found. Subsequently, the extracted pesticide-hThr83, pesticide-hTyr124, pesticide-hTyr341, and pesticide-hHis447 HB forming geometries were QM optimized in order to locate the transition states (TSs) and/or the intermediary states (ISs) by which the protons transfers occur in the particular segment of biochemical reaction. The second level of calculation attempted to find the identical TSs and/or ISs upon the manual adjustment of starting pesticide-hThr83, pesticide-hTyr124, pesticide-hTyr341, and pesticide-hHis447 geometries, by means of the TS-generating hydrogen atom position between the pesticide and the regarded AChE residue, in a kind of validation process. Each of the reaction pathways assumed the contemplation by means of Intrinsic Reaction Coordinate (IRC) analysis. Remarkably, the independent calculations gave similar/identical results, confirming that each of the obtained TSs and/or ISs is not formed by chance.
Thus, the optimization of atrazine-hThr83, atrazine-hTyr124, atrazine-hTyr341, and atrazine-hHis447 HB geometries confirmed that the proton transfer, within each of the geometries, occurs via the adequate transition states. Hence, the inhibition of hAChE starts with the nucleophilic attack of atrazine’s m-N′ atom hydrogen towards hThr83 (Figure 8a). The distance between the hThr83 hydroxyl group oxygen atom as a HB acceptor and m-N′ atom as a HB donor amounts 3.758 Å, characterizing particular HB as weak, almost electrostatic by character. This is, certainly, the expected increase in atrazine-hThr83 HB distance, arisen following the QM optimization, in comparison to the one recorded during the MD simulations (dHB = 2.253−2.976 Å). Nevertheless, even the weak HB is enough to increase the acute toxicity of atrazine to humans, in comparison to the propazine and simazine, inasmuch as the latter pesticides are, according to the QM studies, unable to donate proton for the formation of described TS1. Instead, the propazine-hThr83 and simazine-hThr83 interactions (Figure S43a and Figure S44a, respectively) were characterized by the IS1; the optimized QM distance between the hThr83 hydroxyl group oxygen atom and propazine and simazine m-N′ atom was 4.259 Å and 4.281 Å, respectively, suggesting the electrostatic nature of the established bond instead of the HB one. Back to the atrazine, upon the formation of atrazine-hThr83 TS1 (Figure 8a and Scheme 1), the deprotonation of hThr83 side chain hydroxyl group occurs, yielding the negatively charged oxygen atom, whereas at the same time m-N′ atom becomes positively charged. Identical TS1 conformation and reaction pathway was obtained upon the manual setup and the optimization of the transition state geometry. The calculated activation energy barrier for the generated TS1 is 11.2 kcal/mol and this initial reaction is the rate-limiting step for the hAChE inhibition by atrazine (Figure 8a and Figure 9). The Free energy profile of all of the described reaction pathways was presented similarly to the study describing the catalytic reaction mechanism of AChE on the ACh level [82]. The manual setup of TS1 coordinates and its optimization resulted in the fact that the 47% of IRC described the proton transfer from TS1 to hThr83, whereas the 53% or IRC quantified the proton transfer from TS1 to atrazine m-N′ atom. This result further supported the proposed nucleophilic attack; the proton transfer is concerted, and the characterized transition state is meaningful. The precise atrazine-hThr83 TS1/hydrogen bond is formed and is stable at all transition and intermediate states during the inhibition process.
The atrazine-based hAChE inhibition further flows towards the deprotonation of m-N′ atom (Figure 8b and Scheme 1) where hTyr341 phenoxy anion serves as a HB-acceptor. The deprotonation concerted via the TS2, in which the distance between the electronegative atoms amounts 3.211 Å. The similar transition states (this time the TS1s) were also observed for propazine and simazine (Figure S43b and Figure S44b, respectively), and for peculiar pesticides this was the rate-limited reaction (Figure S45a,b). The calculated activation energy barriers for the generation of propazine-hThr83 and simazine-hThr83 TS1 complexes were 8.5 and 8.9 kcal/mol, respectively (Figure S45a and Figure S45b, respectively). To adopt the atrazine-hTyr341 TS2 geometry (Figure 8b and Scheme 1), atrazine is slightly in-plane rotated towards hTyr341, whereas propazine and simazine suffered more severe longitudinal movement (Figure S43b and Figure S44b, respectively). Considering that the energy difference between the atrazine-hTyr341 reactant state and TS2 is only 1.6 kcal/mol (Figure 9), an extra stabilization comes with the formation of atrazine-hTyr341 TS2. The length of atrazine-hTyr124 HB characterizes it as a moderate one. The second level of QM optimization also confirmed the meaningfulness of TS2. As previous, the meaningfulness of TS2 is confirmed since the 56% of IRC described the proton gliding from TS3 to hTyr341, whereas the 44% or IRC illustrated the proton transfer form TS2 for atrazine m-N′ atom. Similar trends were noticed for propazine-hTyr341 and simazine-hTyr341 TS1 (Figure S43b and Figure S44b, respectively). Nevertheless, during the rotation of atrazine (the movement of propazine or simazine), the TS1 (IS1) is sustained (Figure S43b and Figure S44b, respectively).
In the third step of the reaction pathway, the electrophilic nature of the re-established secondary m-N′ amine is expressed by the release of proton towards the hTyr124 hydroxyl group anion. The proton transfer takes place via the TS3, in which the distance between the electronegative atoms amounts 3.522 Å (Figure 8c and Scheme 1). For the purpose of TS3 generation, atrazine is forced to move longitudinally in order to form the tetrahedral structure with the HB-involved residues. Therefore, this step is slightly energetically more expensive, with the free energy difference of 4.8 kcal/mol (Figure 9). The characteristic of the tetrahedral TS3 structure is that it holds the negatively charged secondary m-N′ nitrogen atom. The distinct state is supported by the aromatic nature of the atrazine’s 1,3,5-triazine ring, which can accept the m-N′ amine negative charge by the intramolecular resonance stabilization. As mentioned above, the meaningfulness of the TS3 is confirmed since 57% of IRC described the proton gliding from TS3 to hTyr124, whereas the 43% or the IRC illustrated the proton transfer from TS3 to atrazine’s m-N′ atom. On the other hand, the propazine-hTyr124 and the simazine-hTyr124 interactions are characterized by the IS2 (Figure S43c and Figure S44c, respectively), in which the QM optimized distances between the hTyr124 and the propazine and simazine were 4.882 and 3.862 Å, respectively.
Moreover, the formation of the atrazine-hAChE tetrahedral TS3 structure appears to be the pre-condition for hHis447 to be jointly arranged with the atrazine (Figure 8d and Scheme 1). In this last step of the reaction pathway, the atrazine’s m-N″ amine suffers the nucleophilic attack from hHis447 imidazole ring, in which the proton transfer is carried out from the pesticide towards the amino acid via the TS4 (the HB distance was 2.699 Å). Consequently, there is a localization of the negative charge around the m-N″ amine as well, which can be once again efficiently delocalized by the atrazine′s 1,3,5-triazine ring. The similar pathway is observed for the propazine-hHis447 and simazine-hTyr447 interaction, with the transition state labelled as TS2 (the HB distances were 2.583, and 2.540 Å, respectively). In comparison with the TS4 (TS2) optimized MD snapshot, the orientation of hHis447 needs to be adjusted very slightly, following which there is a proton transfer from the atrazine/propazine/simazine m-N″ amine to the residue accompanied with the spontaneous break of the scissile m-N″-H bond (Figure 8d, Figure S43d and Figure S44d, respectively). As a result, a short and low-barrier hydrogen bond (LBHB) is formed between the hHis447 and the atrazine’s/propazine’s/simazine’s m-N″ amine (Figure 8d, Figure S43d and Figure S44d, respectively), and the willing proton transfer is thus observed. Energetically, this is the most favourable reaction since the free energy differences between the reactant and transition states were 0.70, 0.68, and 0.62 kcal/mol, respectively. Alongside with the TS1, the TS4 (TS2) is by far the most important transition state as it interrupts the hHis447 to serve as a base for hAChE catalytic triad. By all means, the meaningfulness of TS4 is confirmed since the 51% of IRC described the proton gliding from TS3 for hHis447, whereas the 49% or IRC illustrated the proton transfer form TS4 for atrazine m-N′ atom (Figure 8d and Scheme 1). Similar trends were also noticed for propazine-hTyr341 and simazine-hTyr341 TS2 (Figure S43d and Figure S44d, respectively). Nevertheless, during the formation of the atrazine-hHis447 TS4 (propazine-hHis447 and simazine-hTyr447 TS2), the TS1 (IS1), TS2 (TS1), and TS3 (IS2) were sustained.
The quantum mechanical studies do not provide the basis for the hSer203 covalent modification by atrazine’s, propazine’s or simazine’s m-N″ atom substituent, despite the fact that in metabolic pathways pesticides may release the particular functional group. The literature survey supports particular findings since the metabolites from atrazine, simazine and propazines are not able to establish the interactions with acetylcholinesterase [83]. This statement was further confirmed in the concentration dependent kinetic studies.

2.8. The Concentration-Dependent Kinetics of Human Acetylcholinesterase Inhibition by 2-Chloro-1,3,5-triazine-based Pesticides

There are numerous evidences that organophosphates are irreversible AChE inhibitors whereas the carbamates acts as reversible AChE inhibitors [73]. However, studies show that atrazine decreases the AChE catalytic activity in Chironomus tentans [84] and Carassius auratus, applied synergistically with organophosphate insecticides, [80]. Nevertheless, the data describing the atrazine’s and the related compounds individually influence the hAChE catalytic activity is limited. Therefore, the 2-chloro-1,3,5-triazine-based pesticides are herein evaluated as hAChE inhibitors, in the concentrations comparable to the calculated LD50 values for Homo sapiens (Figure 10).
The results show that all of the examined 2-chloro-1,3,5-triazine-based pesticides inhibit Homo sapiens AChE activity in a concentration-dependent manner (Figure 10). The inhibition curves analysis implies that the aging reaction, in which the hAChE releases itself from covalently bound pesticide, occurs in the presence of atrazine, propazine, and simazine at a high rate, and that the regarded pesticides most likely do no inhibit the hAChE covalently [85]. Therefore, the protocol by which the Homo sapiens LD50 values were calculated, as well as all the performed molecular modelling studies of 2-chloro-1,3,5-triazine-based pesticides, were experimentally validated.

3. Materials and Methods

3.1. Materials

Acetylcholine chloride (CAS: 60-31-1), acetylcholinesterase from human erythrocytes (CAS: 9000-81-1), and DTNB (5,5′-dithiobis(2-nitrobenzoic acid)) (CAS: 69-78-3), atrazine (CAS: 1912-24-9), propazine (CAS: 139-40-2), simazine (CAS: 122-34-9), and PBS (MDL: MFCD00131855) were purchased from Sigma Aldrich (St. Louis, Mo., USA).

3.2. The Generation of Pesticides Structures

The crystal structures of simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] were retrieved from The Cambridge Crystallographic Data Centre (CCDC) (CSD IDs: KUYXIM, ULEJIF, IPCPYB, HANBAA). The remaining pesticides structures were generated by drawing using the MacroModel version 8 software (Schrodinger, Cambridge, MA, USA) [75].

3.3. The Conformational Analysis

The conformational analyses were performed by using the MacroModel version 8 software [75] via the Monte Carlo/multiple minimum (MC/MM) approach [76]. The solute was described using five different potentials (MM3, AMBER94, MMFF, MMFFs, OPLSAA) [77,85,86,87,88]. The effect of solvent was incorporated into the MC/MM calculations using the generalized Born/surface area (GB/SA) continuum solvent model [89]. Cut-offs of 12.0 Å and 7.0 Å were employed for electrostatic and van der Waals non-bonded interactions, respectively [90]. The MC simulation involved 104 steps at 310.15 K, applied to all the rotatable bonds, with random torsional rotations of up to ±180°. This was combined with 103 steps of energy minimization. All the conformational calculations were performed using the MacroModel 8.0 and BatchMin suite of programs [91].

3.4. The Mus musculus and Homo sapiens AChE Sequences Alignment

The alignment of mAChE and hAChE sequences retrieved from the UniProt database (entries P221836 and P22303, respectively, Figure S1) was performed by means of ClustalW module [92] using the standard setup.

3.5. The Mus musculus and Homo sapiens AChE Complexes Preparation

To the best of our knowledge, there are no crystal structures containing co-crystallized structures of pesticides listed in Table 1, in complex with either Mus musculus or Homo sapiens AChE. Therefore, for the purpose of targeted pesticides binding mode definition into the active sites of Mus musculus and Homo sapiens AChEs, respectively, 59 acetylcholine esterase-inhibitor complexes (51 containing mice inhibitors and 8 saturated with human inhibitors) were collected from Protein Data Bank and submitted to structure-based alignment assessment protocols (see Tables S9 and S10). Initially, complexes were prepared for molecular modelling after being loaded into the UCSF Chimera v1.10.1 software [93] for Linux 64 bit architecture and visually inspected. The downloaded complexes were homodimers. For the experimental purposes preparation, complexes were initially arbitrarily superimposed using human crystal under the code name 4EY5 [94] as a template (the best-resolved complex with the resolution of 2.3 Å) using the Matchmaker module and then separated in chains using the command line implementation of the Chimera split command. Upon the thorough examination, the B chain of each complex was retained for further manipulation. Compared to chains A containing complexes, those with chains B were more complete by means of inhibitors presence. Inhibitors were extracted from each chain B complexes after which either proteins or ligands were improved by assigning the hydrogen atoms at pH of 7.4 and Amber parameters [95] using Antechamber semi-empirical QM method. Complexes were the regenerated and energy minimized as follows: through the leap module they were solvated with water molecules (TIP3P model, SOLVATEOCT Leap command) in a box extending 10 Å in all directions, neutralized with either Na+ or Cl ions, and refined by a single point minimization using the Sander module with maximum 1000 steps of the steepest-descent energy minimization followed by maximum 4000 steps of conjugate-gradient energy minimization, with a non-bonded cutoff of 5 Å. Minimized complexes were re-aligned (4EY5 as a template), after which all of the ligands were extracted to compose the SB aligned training set ready to be utilized for the subsequent structure-based alignment assessment.

3.6. The Structure-Based Alignment Assessment

The SB alignment assessment for either Mus musculus or Homo sapiens AChE inhibitors was performed by means of the four step standard protocol [78,96], by means of either AutoDock [79], AutoDock Vina [80], and DOCK [81] docking algorithms/scoring functions.
Experimental conformation re-docking (ECRD): a procedure in which the experimental conformations (EC) are flexibly docked back into the corresponding protein, evaluating the program for its ability to reproduce the observed bound conformations.
Randomized conformation re-docking (RCRD): similar assessment to ECRD with a difference that the active site of the protein is virtually occupied by conformations initially obtained from computational random optimization of corresponding co-crystallized molecules coordinates and positions. Here the programs are evaluated for their ability to find the experimental pose starting by a randomized minimized conformation.
Experimental conformation cross-docking (ECCD): comparable to ECRD, but the molecular docking was performed on all the training set proteins except for the corresponding natives. Here the programs are evaluated to find the binding mode of ligand in a similar but different from the native complex active site at the same time mimicking discrete protein flexibility.
Randomized conformation cross-docking (RCCD): same as the ECCD but using RCs as starting docking conformations. This is the highest level of difficulty since the program is demanded to dock any given molecule into an ensemble of protein conformations not containing the native one. The outcome is considered as the most important ability of docking program as the most accurate scoring function is applied to any test set molecules of which the experimental binding mode is unknown. The related docking accuracy (DA) [97] is a direct function of the program’s probability to find a correct binding mode for an active molecule.

3.6.1. The AutoDock Settings

The prepared protein structures were imported into AutoDockTools graphical user interface. Nonpolar hydrogen atoms were removed while Gasteiger charges and solvent parameters were added. All of the tested compounds were used as ligands for docking. The rigid root and rotatable bonds were defined using AutoDockTools. The docking was performed with AutoDock 4.2. The xyz coordinates (in Ångströms) of the cuboid grid box used for the computation were: for mouse inhibitors Xmin/Xmax = 1.420/31.800, Ymin/Ymax = 108.000/152.500, Zmin/Zmax = 8.200/43.200; for human inhibitors Xmin/Xmax = 3.200/29.500, Ymin/Ymax = 109.200/153.800, Zmin/Zmax = 8.800/42.100, with a purpose to embrace all the minimized inhibitors spanning 10 Å in all three dimensions. The Lamarckian Genetic Algorithm was used to generate orientations or conformations of ligands within the binding site. The procedure of the global optimization started with a sample of 200 randomly positioned individuals, a maximum of 1.0 × 106 energy evaluations and a maximum of 27,000 generations. A total of 100 runs were performed with RMS Cluster Tolerance of 0.5 Å.

3.6.2. The Autodock Vina Settings

The docking simulations were carried out in the same grid space with an energy range of 10 kcal/mol and exhaustiveness of 100 with RMS Cluster Tolerance of 0.5 Å. The output comprised 20 different conformations for every receptor considered.

3.6.3. The DOCK6 Settings

During the docking simulations with DOCK program, the proteins were rigid while the inhibitors flexible and subjected to an energy minimization process. The solvent accessible surface of each enzyme without hydrogen atoms was calculated using the DMS program [98], using a probe radius of 1.4 Å. The orientation of ligands was described using the SPHGEN module and by sphere selector. A box around the binding site is constructed with the accessory program SHOWBOX. The steric and electrostatic environment of the pocket was evaluated with the program Grid using a 0.3 Å of grid spacing. Selected spheres were within 8 Å from ligand heavy atoms of the crystal structure and for computing the energy grids an 8 Å box margin and 6–9 VDW exponents were used.

3.7. The Ligand-Based and Structure-Based Alignment Accuracy

The alignment fitness for LB or SB approach was quantified by evaluating the RMSD values. The alignment accuracy was initially quantified after reproducing the known crystal structures of simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12], with the available force fields, to obtain the conformational alignment accuracy (CAA). Following, the superposition of pesticides structures obtained after the conformational analysis was performed with the purpose to compare the force fields performances, i.e., the force filed-dependent alignment accuracy, FFDAA. Finally, superimposition of docking poses provided the docking accuracy, DA. CAA, FFDAA or DA can be used to test how the conformational alignment can predict the ligand pose as close as possible to the experimentally observed one [97]. The alignment accuracy values can be calculated by the following Equation (7):
x A = f r m s d ,   a + 0.5 ( f r m s d   b f r m s d   a )
In particular: xA = CAA or FFDAA, in the case of LB alignment accuracy, whereas the xA = DA, in the case of SB accuracy. The frmsda and frmsdb are fractions of aligned ligands showing an RMSD value less than or equal to 2 Å (a coefficient) and 3 Å (b coefficient), respectively. The widely accepted standard is that the correctly aligned conformations are those with RMSD less than 2 Å on all the heavy atoms between the generated and crystallographic structure, when considering the CAA and FFDAA, while the correctly docked structures are those with RMSD less than 2 Å on all the heavy atoms between docked and co-crystallized ligand, when considering the DA. Structures with the RMSD values larger than 4 Å were considered incorrectly aligned/docked. Structures with RMSD between 2 and 3 were considered partially aligned/docked, whereas those with RMSD higher than 3 and were misaligned/mis-docked and thus not considered in the FFDAA/DA calculation.

3.8. The Structure-Based Alignment of Pesticides

The molecular docking of pesticides within either mAChE and hAChE was conducted by the AutoDock Vina [80]. The molecular docking studies were carried out in the cuboid grid box with the following xyz coordinates (in Ångströms): for mAChE Xmin/Xmax = 1.420/31.800, Ymin/Ymax = 108.000/152.500, Zmin/Zmax = 8.200/43.200; hAChE Xmin/Xmax = 3.200/29.500, Ymin/Ymax = 109.200/153.800, Zmin/Zmax = 8.800/42.100, with the energy range of 10 kcal/mol and exhaustiveness of 100 with RMS Cluster Tolerance of 0.5 Å. The output comprised 20 different conformations for every receptor considered.

3.9. The Pesticide-AChE Complexes Molecular Dynamics Simulations

The pesticide-AChE complexes formed after molecular docking were used to perform the explicit solvent molecular dynamics (MD) simulations. The parallelized Desmond Molecular Dynamics System [99] and the associated analysis tools, available within the Schrödinger suite 2009 [75] were used for this purpose. The models were prepared using the ‘system builder’ utility. The MMFF force field parameters were assigned to all the simulation systems. Each inhibitors-enzyme complex was placed in the centre of an orthorhombic box ensuring 10 Å solvent buffers between the solute and the boundary of the simulation box in each direction. The TIP3P model was used to describe the water molecules. Moreover, Na+ and Cl ions were added at a physiological concentration of 0.15 M. The model systems were relaxed to 0.5 Å using the single step minimization protocol and were subjected to the production phase, with the NPT ensemble and the periodic boundary conditions for 1.2 ns. The pressure control was used to maintain the pressure of the cell (1.013 bars) with the isotropic coupling method. The Nose-Hoover thermostat method was used to control the temperature at 310.15 K. The Heavy atom-hydrogen covalent bonds were constrained using the SHAKE algorithm which allowed 2 fs integration step to be used. The integration of the motion equation using the RESPA multiple time step scheme was used in the calculation. The Long-range electrostatic forces were obtained using the smooth Particle-Mesh Ewald (PME) method. In order to calculate the short-range electrostatics and the Lennard-Jones interaction, the cut-off distance was set to 9.0 Å, and the trajectories and the energies were recorded at every 4.8 ps. The simulation quality analysis tool was used to analyse the trajectories; RMSD and RMSF values, the hydrogen bond distances, the angles, and the van der Waals interactions were obtained over the simulation trajectories.

3.10. The MM-GBSA Calculations and Free Energy Decomposition

The binding free energy [100] of each system was calculated using the MM-GBSA calculations according to the following Equation (8):
                  Δ G bind   = Δ E ele + Δ E vdW + Δ G solv + T Δ S ,  
where ΔEele corresponds to the electrostatic energy difference between the receptor-ligand bound and the calculated unbound states using the OPLS_2005 force field, ΔEvdW corresponds to the van der Waals contribution, while ΔGsolv is the corresponding solvation free energy contribution of the binding which was calculated using the GB/SA continuum model. The Embrace package implemented in MacroModel was used for the ΔEele, ΔEvdW, and ΔGsolv calculations. The entropy change ΔS was calculated using the Rigid Rotor Harmonic Oscillator (RRHO) calculations. Having used this algorithm, the change in vibrational, rotational, and translational (VRT) entropy of the ligands on binding was estimated. For the RRHO calculations, the representative complexes were pre-minimized using the Desmond with the explicit solvent retained; a 2000 steps LBFGS minimization (first 10 steps steepest descent) with the residues beyond 15 Å of ligands restrained and a convergence criterion of 0.05 kcal mol−1 Å−1 was used.

3.11. The Quantum Mechanical Mechanistic Studies

The DFT-based calculations were performed on the MD-generated pesticide-hAChE complexes on three different levels. The first level included the extraction of hydrogen bond-forming geometries of MD equilibrated pesticide-hThr83, pesticide-hTyr124, pesticide-hTyr341, and pesticide-hHis447 sub-complexes; subsequently, the extracted geometries were QM optimized to locate the transition states (TSs) and/or intermediary states (ISs), by which the protons transfers occur on a particular segment of a biochemical reaction. The second level of the calculation attempted to find the identical TSs and/or ISs following the manual adjustment of the TS-generating hydrogen atom position, between the pesticide and the regarded AChE residue, within the pesticide-hThr83, pesticide-hTyr124, pesticide-hTyr341, and pesticide-hHis447 geometries. Each of the reaction pathways assumed the contemplation by means of Intrinsic Reaction Coordinate (IRC) analysis. As water molecules are not involved in HBs generation, they were all discarded from the system. The initiation or the manual setup of geometries was performed by means of GaussView6 package [101] following which the B3LYP method, as implemented in Gaussian 09 (Gaussian Inc., Wallingford, CT, USA) [102] was used to explore the geometry of the reactants, transition states, intermediary states and products. The single point calculations were carried out on all the geometries at the B3LYP/6-31G* level of theory [82]. The saddle points and transition states were quantified by means of frequency calculations.

3.12. The Assay of the AChE Activity

The AChE catalytic activity was measured at 22 °C by the Ellman method [103]. The assay mixture contained the 0.25 mM ACh, 0.25 mM DTNB, and 50 mM sodium phosphate buffer. The remaining assay conditions have been reported previously. For the selection of a suitable concentration of AChE, at which the relationship between the initial velocity and the total enzyme concentration should be linear, 0.20–3.20 mg/mL of AChE was assayed in vitro from 2 to 20 min under the same conditions of the temperature and the pH. The rate of change of the substrate cleavage, determined by measuring the absorbance of the reaction product at a wavelength of 412 nm, was performed at different time intervals. The product was calculated for each [AChE] at 4-min incubation time intervals. To study the effect of atrazine, simazine, or propazine on AChE substrate cleavage, the enzyme was preincubated with each pesticide at different concentrations ranging 0–6 μg/mL (concentrations corresponding to the LD50 values) for 10 min prior to the addition of the substrate.

4. Conclusions

The present report describes the application of the molecular modelling techniques to shed light on the pharmacology of the commonly used pesticides atrazine, simazine, propazine, carbofuran, monocrotophos, dimethoate, carbaryl, tebufenozide, imidacloprid, acetamiprid, diuron, monuron, and linuron as AChE inhibitors The pesticides’ putative pre-binding conformations in the inter-synaptic environment were predicted by means of the conformational analysis where the determined global minima values and the ability to reproduce the pesticides experimental conformations indicated MMFF as the best performing FF. Having known the pesticides pre-bound conformations, their acute toxicities were interrelated to the Eglob_min values through the QSAR studies: for the majority of pesticides, a clear inverse correlation could be observed between the Eglob_min values and the high acute toxicities. However, the chemometric analysis indicated that the Eglob_min values could not be used as proper indicators for high or low LD50 values against Mus musculus or Homo sapiens.
Therefore, the pesticides acute toxicity was further on regarded trough the pesticides bound conformations, upon the molecular docking and molecular dynamics-driven interactions with either Mus musculus or Homo sapiens AChE. The conducted studies confirmed the binding modes and the pharmacology of evaluated pesticides as reversible and reversible-like (i.e., future covalent) AChE inhibitors and the determined pesticides pharmacodynamics origin of acute toxicity. The superposed chemometric investigation revealed that for all the pesticides, a clear interrelationship could be observed between the ∆Gbinding values and high acute toxicity, indicating that the ∆Gbinding values can be used as indexes for high or low LD50 values against Mus musculus. Moreover, the obtained chemometric guideline, in the form of the QSAR model pLD50 = −0.7885∆Gbinding − 2.44, is, by all means, suitable for the prognosis of LD50 values against Homo sapiens and can be further used as a universal tool for the acute toxicities administration of novel pesticides. Further studies which consider the inclusion of additional parameters and machine learning algorithm application are in due course in order to generate quantitative structure-activity relationships models able to include all considered pesticides in a comprehensive mathematical model.
Nevertheless, it is interesting that the global minima or free energy of formation values may indicate the level of acute toxicity; the calculations may discard the potential new pesticides even before the actual application in agriculture. Therefore, the correlations between the pLD50 and Eglob_min or ∆Gbinding are hereby presented for the very first time as an unprecedented way to study pesticides and to predict their acute toxicity.
Moreover, the subsequent structure-based quantum mechanical mechanistic studies for 2-chloro-1,3,5-triazine-based pesticides, as the least understood group by means of pharmacology, revealed pathways by which considered pesticides poison Homo sapiens AChE in a reversible fashion manner. To summarize, all the noticed interactions could further be used in the discovery of novel pesticides with the desirable lower acute toxicity against humans. It is important to emphasize that the conducted QM studies were performed to provide the structural basis for the pesticides’ human AChE toxicity, not the general toxicity. The pesticide toxicity is mainly carried out by means of mechanism; general toxicity is normally due to off target interactions or to active metabolites.

Supplementary Materials

The Supplementary Materials are available online.

Author Contributions

B.B.A. and M.M. designed the study; B.B.A., A.R., J.S.M., T.M.T.-P. and R.M. performed conformational analyses. M.M., N.S., N.M., and R.R. performed statistical studies, molecular docking, molecular dynamics studies and QM DFT studies.

Funding

This research was funded by the Ministry of Education and Science of Republic of Serbia (Grant Nos. III45006, 174007, and III43004) and Progetti di Ricerca di Università 2015, Sapienza Università di Roma (Grant Nos. C26A15RT82 and C26A15J3BB).

Acknowledgments

Gratitude is expressed to Goran Bogdanovic, Institute of Nuclear Sciences Vinča, Belgrade, Republic of Serbia, for providing pesticides crystal structures. Special thanks are due to Jill Barber, Division of Pharmacy and Optometry, University of Manchester, Manchester, UK for providing computational facilities and all comments and suggestions which improved the quality of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Gilden, R.C.; Huffling, K.; Sattler, B. Pesticides and health risks. J. Obstet. Gynecol. Neonatal Nurs. 2010. [Google Scholar] [CrossRef] [PubMed]
  2. Arduini, F.; Ricci, F.; Tuta, C.S.; Moscone, D.; Amine, A.; Palleschi, G. Detection of carbamic and organophosphorous pesticides in water samples using a cholinesterase biosensor based on Prussian Blue-modified screen-printed electrode. Anal. Chim. Acta 2006. [Google Scholar] [CrossRef] [PubMed]
  3. U.S. Environmental Protection Agency, Washington, DC (1994) Pesticide Industry Sales and Usage, 1992 and 1993 Market Estimates. Available online: https://nepis.epa.gov/Exe/ZyPURL.cgi?Dockey=200001EL.TXT (accessed on 15 June 2018).
  4. Bolton, T.B.; Lim, S.P. Action of acetylcholine on smooth muscle. Z. Kardiol. 1991, 80 (Suppl. 7), 73–77. [Google Scholar] [PubMed]
  5. Jones, B.E. From waking to sleeping: Neuronal and chemical substrates. Trends Pharmacol. Sci. 2005, 26, 578–586. [Google Scholar] [CrossRef] [PubMed]
  6. Chapalamadugu, S.; Chaudhry, G.R. Microbiological and biotechnological aspects of metabolism of carbamates and organophosphates. Crit. Rev. Biotechnol. 1992. [Google Scholar] [CrossRef] [PubMed]
  7. Đorđević, J.; Vladisavljević, G.T.; Trtić-Petrović, T. Liquid-phase membrane extraction of targeted pesticides from manufacturing wastewaters in a hollow fibre contactor with feed-stream recycle. Environ. Technol. 2017, 38, 78–84. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  8. Đorđević, J.S.; Vladisavljević, G.T.; Trtić-Petrović, T.M. Removal of the selected pesticides from a water solution by applying hollow fiber liquid–liquid membrane extraction. Ind. Eng. Chem. Res. 2014, 53, 4861–4870. [Google Scholar] [CrossRef] [Green Version]
  9. Brown, J.E.; Baughman, R.G. N-Ethyl-6-ethylamino-4-oxo-1,3,5-triazin-2-aminium chloride (Oxysimazine·HCl). Acta Crystallogr. Sect. E 2010. [Google Scholar] [CrossRef] [PubMed]
  10. Cheon, S.; Kim, T.H.; Park, K.-M.; Kim, J. Monocrotophos: Dimethyl (E)-1-methyl-2-(methylcarbamoyl)ethenyl phosphate. Acta Crystallogr. Sect. E 2011, 67, o584. [Google Scholar] [CrossRef] [PubMed]
  11. Lapp, R.L.; Jacobson, R.A. Crystal and molecular structures of organophosphorus insecticides. 13. S-isopropyl O-methyl O-(3,5,6-trichloro-2-pyridyl) phosphoramidothioate and dimethoate. J. Agric. Food Chem. 1980, 28, 755–759. [Google Scholar] [CrossRef]
  12. Chopra, D.; Mohan, T.P.; Rao, K.S.; Guru Row, T.N. IUCr (E)-N1-[(6-Chloropyridin-3-yl)methyl]-N2-cyano-N1-methylacetamidine. Acta Crystallogr. Sect. E 2004, 60, o2374–o2375. [Google Scholar] [CrossRef]
  13. Boudina, A.; Emmelin, C.; Baaliouamer, A.; Grenier-Loustalot, M.F.T.; Chovelon, J.M. Photochemical behaviour of carbendazim in aqueous solution. Chemosphere 2003. [Google Scholar] [CrossRef]
  14. World Health Organization. The WHO Recommended Classification of Pesticides by Hazard and Guidelines to Classification 2009; World Health Organization: Geneva, Switzerland, 2010; pp. 1–60. [Google Scholar]
  15. Bravo, R.; Caltabiano, L.M.; Weerasekera, G.; Whitehead, R.D.; Fernandez, C.; Needham, L.L.; Bradman, A.; Barr, D.B. Measurement of dialkyl phosphate metabolites of organophosphorus pesticides in human urine using lyophilization with gas chromatography-tandem mass spectrometry and isotope dilution quantification. J. Expo. Anal. Environ. Epidemiol. 2004. [Google Scholar] [CrossRef] [PubMed]
  16. Catalá-Icardo, M.; López-Paz, J.L.; Choves-Barón, C.; Peña-Bádena, A. Native vs. photoinduced chemiluminescence in dimethoate determination. Anal. Chim. Acta 2012. [Google Scholar] [CrossRef] [PubMed]
  17. Park, C. A Dictionary of Environment and Conservation, 1st ed.; Oxford University Press: Oxford, UK, 2007. [Google Scholar]
  18. Elbert, A.; Haas, M.; Springer, B.; Thielert, W.; Nauen, R. Applied aspects of neonicotinoid uses in crop protection. Pest Manag. Sci. 2008, 64, 1099–1105. [Google Scholar] [CrossRef] [PubMed]
  19. USEPA. Reregistration Eligibility Decision (RED), Linuron. 1995. Available online: https://archive.epa.gov/pesticides/reregistration/web/pdf/0047.pdf (accessed on 15 June 2018).
  20. Gallo, M.A.; Lawryk, N.J. Organic phosphorus pesticides. In Handbook of Pesticide Toxicology, Hayes, W.J., Laws, E.R.J., Eds.; Academic Press: San Diego, CA, USA, 1991. [Google Scholar]
  21. U.S. Environmental Protection Agency. Propazine: Health Advisory. Office of Drinking Water, US EPA, Washington, DC, USA. 1988. Available online: https://www.epa.gov (accessed on 16 March 2016).
  22. Extension Toxicology Network. Pesticide Information Profiles, Oregon State University, Diuron. 1996. Available online: http:// www.extoxnet.orst.edu (accessed on 19 March 2016).
  23. Worthing, C.R.; Walker, S.B. The Pesticide Manual: A World Compendium, 8th ed.; The British Crop. Protection Council: Thornton Heath, UK, 1987. [Google Scholar]
  24. U.S. Environmental Protection Agency. Simazine: Health advisory. Office of Drinking Water, US EPA, Washington, DC, USA. August 1987. Available online: https://www.epa.gov (accessed on 16 March 2016).
  25. Hayes, W.J.J.; Laws, E.R.J. Handbook of Pesticide Toxicology; Academic Press: New York, NY, USA, 1991. [Google Scholar]
  26. Rotterdam Convention—Operation of the Prior Informed Consent (PIC) Procedure for Banned or Severely Restricted Chemicals in International Trade, Decision Guidance Document, Monocrotophos. Available online: http://www.pic.int (accessed on 16 March 2016).
  27. Occupational Health Services, Inc. MSDS for Dimethoate. OHS Inc., Secaucus, NJ, USA. 1991. Available online: https://hr.umich.edu (accessed on 17 March 2016).
  28. U.S. Environmental Protection Agency, Office of Drinking Water (1987) Carbaryl Health Advisory. Draft Report. Available online: https://www.epa.gov (accessed on 16 March 2016).
  29. Morrison, R.D.; Hamilton, J.D. RH-75992 Technical: Acute Oral Toxicity Study in Male and Female Mice. Unpublished Report No. 91R-003. Rohm Haas Co., Toxicology Department: Philadelphia, PA, USA. Submitted to WHO by Rohm & Haas Co., Philadelphia, PA, USA. Available online: www.inchem.org/documents/jmpr/jmpmono/v96pr11.htm (accessed on 14 June 2018).
  30. Kidd, H.; James, D. Agrochemicals Handbook, 3rd ed.; Royal Society of Chemistry: Cambridge, UK, 1991. [Google Scholar]
  31. O’Neil, M.J. The Merck Index—An Encyclopedia of Chemicals, Drugs, and Biologicals, 13th ed.; Merck Research Laboratories: Whitehouse Station, NJ, USA, 2001; Available online: https://www.rsc.org/merck-index (accessed on 19 March 2016).
  32. Reagan-Shaw, S.; Nihal, M.; Ahmad, N. Dose translation from animal to human studies revisited. FASEB J. 2008, 22, 659–661. [Google Scholar] [CrossRef] [PubMed]
  33. Davies, J.E.; Freed, V.H. An Agromedical Approach to Pesticide Management Some Health and Environmental Considerations; Consortium for International Crop. Protection: Berkeley, CA, USA, 1981. [Google Scholar]
  34. Carvalho, A.T.P.; Barrozo, A.; Doron, D.; Vardi Kilshtain, A.; Major, D.T.; Kamerlin, S.C.L. Challenges in computational studies of enzyme structure, function and dynamics. J. Mol. Graph. Model. 2014. [Google Scholar] [CrossRef] [PubMed]
  35. Jensen, J.S.; Jørgensen, F.S.; Klemmensen, P.D.; Hacksel, U.; Pettersson, I. Conformational analysis of the fungicide fenpropimorph by molecular mechanics calculations and NMR spectroscopy. Pestic. Sci. 1992, 36, 309–318. [Google Scholar] [CrossRef]
  36. Oakeshott, J.G.; Devonshire, A.L.; Claudianos, C.; Sutherland, T.D.; Horne, I.; Campbell, P.M.; Ollis, D.L.; Russell, R.J. Comparing the organophosphorus and carbamate insecticide resistance mutations in cholin-and carboxyl-esterases. Chem. Biol. Interact. 2005, 157, 269–275. [Google Scholar] [CrossRef] [PubMed]
  37. Rathod, A.L.; Garg, R.K. Chlorpyrifos poisoning and its implications in human fatal cases: A forensic perspective with reference to Indian scenario. J. Forensic Leg. Med. 2017. [Google Scholar] [CrossRef] [PubMed]
  38. Sabater, C.; Carrasco, J.M. Effects of the organophosphorous insecticide fenotrothion on growth in five freshwater species of phytoplankton. Environ. Toxicol. 2001, 16, 314–320. [Google Scholar] [CrossRef] [PubMed]
  39. Buratti, F.M.; Volpe, M.T.; Meneguz, A.; Vittozzi, L.; Testai, E. CYP-specific bioactivation of four organophosphorothioate pesticides by human liver microsomes. Toxicol. Appl. Pharmacol. 2003, 186, 143–154. [Google Scholar] [CrossRef]
  40. Howard, J.J.; Oliver, J. Impact of naled (Dibrom 14) on the mosquito vectors of eastern equine encephalitis virus. J. Am. Mosq. Control Assoc. 1997, 13, 315–325. [Google Scholar] [PubMed]
  41. Davis, M.K.; Boone, J.S.; Moran, J.E.; Tyler, J.W.; Chambers, J.E. Assessing intermittent pesticide exposure from flea control collars containing the organophosphorus insecticide tetrachlorvinphos. J. Expo. Sci. Environ. Epidemiol. 2008, 18, 564–570. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Garcia, S.; Abu-Qare, A.; Meeker-O’Connell, W.; Borton, A.; Abou-Donia, M. Methyl parathion: A review of health effects. J. Toxicol. Environ. Health Part B 2003, 6, 185–210. [Google Scholar] [CrossRef] [PubMed]
  43. Dahlgren, J.G.; Takhar, H.S.; Ruffalo, C.A.; Zwass, M. Health effects of diazinon on a family. J. Toxicol. Clin. Toxicol. 2004, 42, 579–591. [Google Scholar] [CrossRef] [PubMed]
  44. Committee for Veterinary Medicinal Products Azamethiphos Summary Report (2). 1999. Available online: http://www.ema.europa.eu/docs/en_GB/document_library/Maximum_Residue_Limits_-_Report/2009/11/WC500010779.pdf (accessed on 14 June 2018).
  45. Bonner, M.R.; Williams, B.A.; Rusiecki, J.A.; Blair, A.; Beane Freeman, L.E.; Hoppin, J.A.; Dosemeci, M.; Lubin, J.; Sandler, D.P.; Alavanja, M.C.R. Occupational exposure to terbufos and the incidence of cancer in the Agricultural Health Study. Cancer Causes Control 2010, 21, 871–877. [Google Scholar] [CrossRef] [PubMed]
  46. Myers, J.P.; Antoniou, M.N.; Blumberg, B.; Carroll, L.; Colborn, T.; Everett, L.G.; Hansen, M.; Landrigan, P.J.; Lanphear, B.P.; Mesnage, R.; et al. Concerns over use of glyphosate-based herbicides and risks associated with exposures: A consensus statement. Environ. Health 2016, 15, 19. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  47. Francis, S.; Saavedra-Rodriguez, K.; Perera, R.; Paine, M.; Black, W.C.; Delgoda, R.; Delgoda, R. Insecticide resistance to permethrin and malathion and associated mechanisms in Aedes aegypti mosquitoes from St. Andrew Jamaica. PLoS ONE 2017, 12, e0179673. [Google Scholar] [CrossRef]
  48. Trachantong, W.; Saenphet, S.; Saenphet, K.; Chaiyapo, M. Lethal and sublethal effects of a methomyl-based insecticide in Hoplobatrachus rugulosus. J. Toxicol. Pathol. 2017, 30, 15–24. [Google Scholar] [CrossRef] [PubMed]
  49. Hofmeister, M.V.; Bonefeld-Jørgensen, E.C. Effects of the pesticides prochloraz and methiocarb on human estrogen receptor α and β mRNA levels analyzed by on-line RT-PCR. Toxicol. In Vitro 2004, 18, 427–433. [Google Scholar] [CrossRef] [PubMed]
  50. Kennedy, G.L. Acute toxicity studies with oxamyl. Fundam. Appl. Toxicol. 1986, 6, 423–429. [Google Scholar] [CrossRef]
  51. Turusov, V.; Rakitsky, V.; Tomatis, L. Dichlorodiphenyltrichloroethane (DDT): Ubiquity, persistence, and risks. Environ. Health Perspect. 2002, 110, 125–128. [Google Scholar] [CrossRef] [PubMed]
  52. Legeay, S.; Clere, N.; Hilairet, G.; Do, Q.-T.; Bernard, P.; Quignard, J.-F.; Apaire-Marchais, V.; Lapied, B.; Faure, S. The insect repellent N,N-diethyl-m-toluamide (DEET) induces angiogenesis via allosteric modulation of the M3 muscarinic receptor in endothelial cells. Sci. Rep. 2016, 6, 28546. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  53. Phosmet. Available online: http://pmep.cce.cornell.edu/profiles/extoxnet/metiram-propoxur/phosmet-ext.html (accessed on 7 April 2018).
  54. Sanderson, D.M. Treatment of poisoning by anticholinesterase insecticides in the rat. J. Pharm. Pharmacol. 1961, 13, 435–442. [Google Scholar] [CrossRef] [PubMed]
  55. Reregistration Eligibility Decision (RED). Available online: https://www3.epa.gov/pesticides/chem_search/reg_actions/reregistration/red_PC-080301_1-Apr-98.pdf (accessed on 12 June 2018).
  56. Chlorpyrifos; No Health Council of The Netherlands: Committee on Updating of Occupational Exposure Limits. Chlorpyrifos; Health-Based Reassessment of Administrative Occupational Exposure Limits. The Hague: Health Council of The Netherlands, 2003; 2000/15OSH/067. 2000, 067; Available online: https://www.gezondheidsraad.nl/sites/default/files/0015067osh.pdf (accessed on 15 June 2018).
  57. Extoxnet Pip Terbufos. Available online: http://extoxnet.orst.edu/pips/terbufos.htm (accessed on 4 July 2018).
  58. Dichlorvos NIOSH Publications and Products. Available online: https://www.cdc.gov/niosh/idlh/62737.html (accessed on 2 July 2018).
  59. Methiocarb Chemical Profile 6/86. Available online: http://pmep.cce.cornell.edu/profiles/insect-mite/fenitrothion-methylpara/methiocarb/methio_prf_0686.html (accessed on 4 July 2018).
  60. Diazinon. Available online: http://www.herbiguide.com.au/Descriptions/hg_DIAZINON.htm (accessed on 2 July 2018).
  61. Ellenhorn, M.J.; Schonwald, S.; Ordog, G.J.; Wasserberger, J. Ellenhorn’s Medical Toxicology: Diagnosis and Treatment of Human Poisoning; Williams Wilkins: Baltimore, MD, USA, 1997; ISBN 068323093X. [Google Scholar]
  62. Pearce, E.M. Kirk-OTHMER encyclopedia of Chemical Technology, 3rd ed.; Wiley-Interscience: New York, NY, USA, 1978. [Google Scholar]
  63. Extoxnet Pip Oxamyl. Available online: http://extoxnet.orst.edu/pips/oxamyl.htm (accessed on 4 July 2018).
  64. “Glyphosate Technical Fact Sheet (Revised June 2015)”. National Pesticide Information Center. 2010. Available online: http://npic.orst.edu/factsheets/archive/glyphotech.html (accessed on 10 May 2018).
  65. Material Safety Data Sheet Malathion sc-211768 Section 1—Chemical Product and Company Identification NFPA Supplier. Available online: https://www.caymanchem.com/msdss/22998m.pdf (accessed on 11 April 2018).
  66. US EPA 2,4-D Reregistration Eligibility Decision, 2005. Associated RED Fact Sheet Archived 2008-05-17 at the Wayback Machine. EPA. Available online: https://archive.epa.gov/pesticides/reregistration/web/html/24d_fs.html (accessed on 14 April 2018).
  67. Methyl Parathion. Available online: http://pmep.cce.cornell.edu/profiles/extoxnet/haloxyfop-methylparathion/methyl-parathion-ext.html#2 (accessed on 2 July 2018).
  68. Dicamba (Banvel) Herbicide Profile 10/83. Available online: http://pmep.cce.cornell.edu/profiles/herb-growthreg/dalapon-ethephon/dicamba/herb-prof-dicamba.html (accessed on 4 July 2018).
  69. CDC-Immediately Dangerous to Life or Health Concentrations (IDLH): Dimethyl-1,2-dibromo-2,2-dichlorethyl Phosphate (Naled)-NIOSH Publications and Products. Available online: https://www.cdc.gov/niosh/idlh/300765.html (accessed on 2 July 2018).
  70. Tsitsanou, K.E.; Thireou, T.; Drakou, C.E.; Koussis, K.; Keramioti, M.V.; Leonidas, D.D.; Eliopoulos, E.; Iatrou, K.; Zographos, S.E. Anopheles gambiae odorant binding protein crystal complex with the synthetic repellent DEET: Implications for structure-based design of novel mosquito repellents. Cell. Mol. Life Sci. 2012, 69, 283–297. [Google Scholar] [CrossRef] [PubMed]
  71. CDC-Immediately Dangerous to Life or Health Concentrations (IDLH): Parathion-NIOSH Publications and Products. Available online: https://www.cdc.gov/niosh/idlh/56382.html (accessed on 2 July 2018).
  72. Department of Pesticide Regulation, C. Summary of Toxicology Data-Sulfoxaflor. 2014. Available online: https://www.cdpr.ca.gov/docs/risk/toxsums/pdfs/6109.pdf (accessed on 5 June 2018).
  73. Čolović, M.B.; Krstić, D.Z.; Lazarević-Pasti, T.D.; Bondzić, AM.; Vasić, V.M. Acetylcholinesterase Inhibitors: Pharmacology and toxicology. Curr. Neuropharmacol. 2013. [Google Scholar] [CrossRef]
  74. Li, Q.; Kong, X.; Xiao, Z.; Zhang, L.; Wang, F.; Zhang, H.; Li, Y.; Wang, Y. Structural determinants of imidacloprid-based nicotinic acetylcholine receptor inhibitors identified using 3D-QSAR, docking and molecular dynamics. J. Mol. Model. 2012, 18, 2279–2289. [Google Scholar] [CrossRef] [PubMed]
  75. Mohamadi, F.; Richards, N.G.J.; Guida, W.C.; Liskamp, R.; Lipton, M.; Caufield, C.; Chang, G.; Hendrickson, T.; Still, W.C. Macromodel: An integrated software system for modeling organic and bioorganic molecules using molecular mechanics. J. Comput. Chem. 1990, 11, 440–467. [Google Scholar] [CrossRef]
  76. Chang, G.; Guida, W.C.; Still, W.C. An internal-coordinate Monte Carlo method for searching conformational space. J. Am. Chem. Soc. 1989. [Google Scholar] [CrossRef]
  77. Allinger, N.L.; Yuh, Y.H.; Lii, J.H. Molecular mechanics. The MM3 force field for hydrocarbons. 1. J. Am. Chem. Soc. 1989, 111, 8551–8566. [Google Scholar] [CrossRef]
  78. Mladenović, M.; Patsilinakos, A.; Pirolli, A.; Sabatino, M.; Ragno, R. Understanding the molecular determinant of reversible human monoamine oxidase B inhibitors containing 2H-chromen-2-one core: Structure-based and ligand-based derived three-dimensional quantitative structure-activity relationships predictive models. J. Chem. Inf. Model. 2017. [Google Scholar] [CrossRef] [PubMed]
  79. Morris, G.M.; Huey, R.; Lindstrom, W.; Sanner, M.F.; Belew, R.K.; Goodsell, D.S.; Olson, A.J. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J. Comput. Chem. 2009, 16, 2785–2791. [Google Scholar] [CrossRef] [PubMed]
  80. Trott, O.; Olson, A. AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization and multithreading. J. Comput. Chem. 2010. [Google Scholar] [CrossRef] [PubMed]
  81. Lang, P.T.; Brozell, S.R.; Mukherjee, S.; Pettersen, E.F.; Meng, E.C.; Thomas, V.; Rizzo, R.C.; Case, D.A.; James, T.L.; Kuntz, I.D. DOCK 6: Combining techniques to model RNA-small molecule complexes. RNA 2009, 15, 1219–1230. [Google Scholar] [CrossRef] [PubMed]
  82. Zhou, Y.; Wang, S.; Zhang, Y. Catalytic reaction mechanism of acetylcholinesterase determined by Born−Oppenheimer ab initio QM/MM molecular dynamics simulations. J. Phys. Chem. B 2010, 114, 8817–8825. [Google Scholar] [CrossRef] [PubMed]
  83. Xing, H.; Wang, J.; Li, J.; Fan, Z.; Wang, M.; Xu, S. Effects of atrazine and chlorpyrifos on acetylcholinesterase and Carboxylesterase in brain and muscle of common carp. Environ. Toxicol. Pharmacol. 2010. [Google Scholar] [CrossRef] [PubMed]
  84. Jin-Clark, Y.; Lydy, M.J.; Zhu, K.Y. Effects of atrazine and cyanazine on chlorpyrifos toxicity in Chironomus tentans (Diptera: Chironomidae). Environ. Toxicol. Chem. 2002. [Google Scholar] [CrossRef]
  85. Rosenfeld, C.A.; Sultatos, L.G. Concentration-dependent kinetics of acetylcholinesterase inhibition by the organophosphate paraoxon. Toxicol. Sci. 2006. [Google Scholar] [CrossRef] [PubMed]
  86. Halgren, T.A. Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94. J. Comput. Chem. 1996. [Google Scholar] [CrossRef]
  87. Wang, J.; Wang, W.; Kollman, P.A.; Case, D.A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graph. Model. 2006. [Google Scholar] [CrossRef] [PubMed]
  88. Cornell, W.D.; Cieplak, P.; Bayly, C.I.; Gould, I.R.; Merz, K.M.; Ferguson, D.M.; Spellmeyer, D.C.; Fox, T.; Caldwell, J.W.; Kollman, P.A. A second generation force field for the simulation of sroteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 1995. [Google Scholar] [CrossRef]
  89. Jorgensen, W.L.; Maxwell, D.S.; Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996. [Google Scholar] [CrossRef]
  90. Arsic, B.; Aguilar, J.A.; Bryce, R.A.; Barber, J. Conformational study of tylosin A in water and full assignments of 1H and 13C spectra of tylosin A in D2O and tylosin B in CDCl3. Magn. Reson. Chem. 2017. [Google Scholar] [CrossRef] [PubMed]
  91. Clark Still, W.; Tempczyk, A.; Hawley, R.C.; Hendrickson, T. Semianalytical treatment of solvation for molecular mechanics and dynamics. J. Am. Chem. Soc. 1990. [Google Scholar] [CrossRef]
  92. Larkin, M.A.; Blackshields, G.; Brown, N.P.; Chenna, R.; McGettigan, P.A.; McWilliam, H.; Valentin, F.; Wallace, I.M.; Wilm, A.; Lopez, R.; et al. Clustal W and clustal X version 2.0. Bioinformatics 2007, 23, 2947–2948. [Google Scholar] [CrossRef] [PubMed]
  93. Pettersen, E.F.; Goddard, T.D.; Huang, C.C.; Couch, G.S.; Greenblatt, D.M.; Meng, E.C.; Ferrin, T.E. UCSF Chimera—A visualization system for exploratory research and analysis. J. Comput. Chem. 2004, 25, 1605–1612. [Google Scholar] [CrossRef] [PubMed]
  94. Cheung, J.; Rudolph, M.J.; Burshteyn, F.; Cassidy, M.S.; Gary, E.N.; Love, J.; Franklin, M.C.; Height, J.J. Structures of human acetylcholinesterase in complex with pharmacologically important ligands. J. Med. Chem. 2012, 55, 10282–10286. [Google Scholar] [CrossRef] [PubMed]
  95. Case, D.A.; Betz, R.M.; Botello-Smith, W.; Cerutti, D.S.; Cheatham, T.E., III; Darden, T.A.; Duke, R.E.; Giese, T.J.; Gohlke, H.; Goetz, A.W.; et al. AMBER 2016; University of California: San Francisco, CA, USA, 2016. [Google Scholar]
  96. Ragno, R.; Ballante, F.; Pirolli, A.; Wickersham, R.B.; Patsilinakos, A.; Hesse, S.; Perspicace, E.; Kirsch, G. Vascular endothelial growth factor receptor-2 (VEGFR-2) inhibitors: Development and validation of predictive 3-D QSAR models through extensive ligand- and structure-based approaches. J. Comput. Aided Mol. Des. 2015. [Google Scholar] [CrossRef] [PubMed]
  97. Bursulaya, B.D.; Totrov, M.; Abagyan, R.; Brooks, C.L. Comparative study of several algorithms for flexible ligand docking. J. Comput. Aided Mol. Des. 2003, 17, 755–763. [Google Scholar] [CrossRef] [PubMed]
  98. Resource for Biocomputing, Visualization, and Informatics (RBVI). Available online: http://rbvi.ucsf.edu (accessed on 6 August 2018).
  99. Bowers, K.J.; Chow, D.E.; Xu, H.; Dror, R.O.; Eastwood, M.P.; Gregersen, B.A.; Klepeis, J.L.; Kolossvary, I.; Moraes, M.A.; et al. Scalable algorithms for molecular dynamics simulations on commodity clusters. In Proceedings of the 2006 ACM/IEEE Conference on Supercomputing, Tampa, FL, USA, 11–17 November 2006; p. 43. [Google Scholar]
  100. Shivakumar, D.; Williams, J.; Wu, Y.; Damm, W.; Shelley, J.; Sherman, W. Prediction of absolute solvation free energies using molecular dynamics free energy perturbation and the opls force field. J. Chem. Theory Comput. 2010. [Google Scholar] [CrossRef] [PubMed]
  101. Dennington, R.; Keith, T.; Millam, J. GaussView 2009. Available online: http://gaussian.com/glossary/g09/ (accessed on 13 March 2018).
  102. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Petersson, G.A.; Nakatsuji, H.; et al. Gaussian 09, Revision A.03; Gaussian, Inc.: Wallingford, CT, USA, 2009. [Google Scholar]
  103. Ellman, G.L.; Courtney, K.D.; Andres, V.; Featherstone, R.M. A new and rapid colorimetric determination of acetylcholinesterase activity. Biochem. Pharmacol. 1961. [Google Scholar] [CrossRef]
Sample Availability: Samples of the compounds are not available from the authors.
Figure 1. The crystal structure of hAChE (PDB ID: 4EY5) in complex with huperzine A (a). Subunit A is depicted in pink ribbons whereas the B subunit is coloured in orange. Huperzine A is presented in black and its structure is highlighted with dashed lines. The anionic sub-site is described by blue spheres, while the esterase sub-area is encircled by red spheres. For the sake of clarity, hydrogen atoms were omitted. The graphic was generated using the UCSF Chimera software (Resource for Biocomputing, Visualization, and Informatics (RBVI), University of California, San Francisco, CA, USA). Biochemical mechanism of ACh hydrolysis by AChE (b).
Figure 1. The crystal structure of hAChE (PDB ID: 4EY5) in complex with huperzine A (a). Subunit A is depicted in pink ribbons whereas the B subunit is coloured in orange. Huperzine A is presented in black and its structure is highlighted with dashed lines. The anionic sub-site is described by blue spheres, while the esterase sub-area is encircled by red spheres. For the sake of clarity, hydrogen atoms were omitted. The graphic was generated using the UCSF Chimera software (Resource for Biocomputing, Visualization, and Informatics (RBVI), University of California, San Francisco, CA, USA). Biochemical mechanism of ACh hydrolysis by AChE (b).
Molecules 23 02192 g001
Figure 2. The plot of experimental vs. fitted pLD50 values calculated with QSAR model 2 for the training set pesticides (training set values are depicted by blue circles; internal test set values are indicated by orange squares) (a); test set pesticides (green circles) (b) against Mus musculus. The plot of calculated vs. fitted pLD50 values calculated with the QSAR model 4 for the training set pesticides (training set values are depicted by blue circles; internal test set values are indicated by orange squares (c); test set pesticides (green circles) (d) against Homo sapiens.
Figure 2. The plot of experimental vs. fitted pLD50 values calculated with QSAR model 2 for the training set pesticides (training set values are depicted by blue circles; internal test set values are indicated by orange squares) (a); test set pesticides (green circles) (b) against Mus musculus. The plot of calculated vs. fitted pLD50 values calculated with the QSAR model 4 for the training set pesticides (training set values are depicted by blue circles; internal test set values are indicated by orange squares (c); test set pesticides (green circles) (d) against Homo sapiens.
Molecules 23 02192 g002
Figure 3. The SB alignment of acetylcholine into mAChE active site (a); hAChE active site (b). The enzyme ribbons are presented in blue and gold, respectively, active site amino acids are depicted in white. For the purpose of clarity, hydrogen atoms are omitted from the presentation. The graphic was generated using the UCSF Chimera software (Resource for Biocomputing, Visualization, and Informatics (RBVI), University of California, San Francisco, CA, USA).
Figure 3. The SB alignment of acetylcholine into mAChE active site (a); hAChE active site (b). The enzyme ribbons are presented in blue and gold, respectively, active site amino acids are depicted in white. For the purpose of clarity, hydrogen atoms are omitted from the presentation. The graphic was generated using the UCSF Chimera software (Resource for Biocomputing, Visualization, and Informatics (RBVI), University of California, San Francisco, CA, USA).
Molecules 23 02192 g003
Figure 4. The SB alignment of atrazine (a); simazine (b); carbofuran (c); monocrotophos (d); dimethoate (e); imidacloprid (f) into the mAChE active site. The enzyme ribbons are presented in blue, the active site amino acids are depicted in white. For the purpose to clarify, the hydrogen atoms are omitted. The graphic was generated using the UCSF Chimera software (Resource for Biocomputing, Visualization, and Informatics (RBVI), University of California, San Francisco, CA, USA).
Figure 4. The SB alignment of atrazine (a); simazine (b); carbofuran (c); monocrotophos (d); dimethoate (e); imidacloprid (f) into the mAChE active site. The enzyme ribbons are presented in blue, the active site amino acids are depicted in white. For the purpose to clarify, the hydrogen atoms are omitted. The graphic was generated using the UCSF Chimera software (Resource for Biocomputing, Visualization, and Informatics (RBVI), University of California, San Francisco, CA, USA).
Molecules 23 02192 g004
Figure 5. The SB alignment of atrazine (a); simazine (b); carbofuran (c); monocrotophos (d); dimethoate (e); imidacloprid (f) into the hAChE active site. The enzyme ribbons are presented in gold; the active site amino acids are depicted in white. For the purpose to clarify, the hydrogen atoms are omitted from the presentation. The graphic was generated using the UCSF Chimera software (Resource for Biocomputing, Visualization, and Informatics (RBVI), University of California, San Francisco, CA, USA).
Figure 5. The SB alignment of atrazine (a); simazine (b); carbofuran (c); monocrotophos (d); dimethoate (e); imidacloprid (f) into the hAChE active site. The enzyme ribbons are presented in gold; the active site amino acids are depicted in white. For the purpose to clarify, the hydrogen atoms are omitted from the presentation. The graphic was generated using the UCSF Chimera software (Resource for Biocomputing, Visualization, and Informatics (RBVI), University of California, San Francisco, CA, USA).
Molecules 23 02192 g005
Figure 6. Schematic view of pesticides binding interactions within the hAChE active site: atrazine, simazine or propazine (a); carbofuran or carbaryl (b); monocrotophos or dimethoate (c) imidacloprid or acetamiprid (d); diuron, monuron, or linuron (e). The hydrogen bonds are presented with pink dashed half arrows, the hydrophobic interactions with a blue dot-line marker, the ionic interactions with the red dashed lines. The proposed mechanism of monocrotophos or dimethoate acute toxicity against hAChE (f).
Figure 6. Schematic view of pesticides binding interactions within the hAChE active site: atrazine, simazine or propazine (a); carbofuran or carbaryl (b); monocrotophos or dimethoate (c) imidacloprid or acetamiprid (d); diuron, monuron, or linuron (e). The hydrogen bonds are presented with pink dashed half arrows, the hydrophobic interactions with a blue dot-line marker, the ionic interactions with the red dashed lines. The proposed mechanism of monocrotophos or dimethoate acute toxicity against hAChE (f).
Molecules 23 02192 g006
Figure 7. The plot of experimental vs. fitted pLD50 values calculated with QSAR model 5 for training set pesticides (a); test set pesticides (b) against Mus musculus. The plot of calculated vs. fitted pLD50 values calculated with QSAR model 5 for training set pesticides (c); test set pesticides (d) against Homo sapiens.
Figure 7. The plot of experimental vs. fitted pLD50 values calculated with QSAR model 5 for training set pesticides (a); test set pesticides (b) against Mus musculus. The plot of calculated vs. fitted pLD50 values calculated with QSAR model 5 for training set pesticides (c); test set pesticides (d) against Homo sapiens.
Molecules 23 02192 g007
Figure 8. The quantum chemical mechanism of Homo sapiens AChE inhibition by atrazine. The extracted geometry of TS1 (a); TS2 (b); TS3 (c); TS4 (d).
Figure 8. The quantum chemical mechanism of Homo sapiens AChE inhibition by atrazine. The extracted geometry of TS1 (a); TS2 (b); TS3 (c); TS4 (d).
Molecules 23 02192 g008
Scheme 1. The mechanism of Homo sapiens AChE inhibition by atrazine, propazine, and simazine, verified by QM studies.
Scheme 1. The mechanism of Homo sapiens AChE inhibition by atrazine, propazine, and simazine, verified by QM studies.
Molecules 23 02192 sch001
Figure 9. Free energy profile for Homo sapiens AChE inhibition by atrazine determined by B3LYP (6-31G*) QM simulations.
Figure 9. Free energy profile for Homo sapiens AChE inhibition by atrazine determined by B3LYP (6-31G*) QM simulations.
Molecules 23 02192 g009
Figure 10. Human AChE inhibition by 2-chloro-1,3,5-triazine-based pesticides.
Figure 10. Human AChE inhibition by 2-chloro-1,3,5-triazine-based pesticides.
Molecules 23 02192 g010
Table 1. Acetylcholine and training set pesticides structures as well as their oral acute toxicities (LD50 values) against Mus musculus (upper values); training set pesticides calculated oral acute toxicities (LD50 values) against Homo sapiens (lower values, bold and italic) [32].
Table 1. Acetylcholine and training set pesticides structures as well as their oral acute toxicities (LD50 values) against Mus musculus (upper values); training set pesticides calculated oral acute toxicities (LD50 values) against Homo sapiens (lower values, bold and italic) [32].
PesticideStructureLD50 (mg/kg)Ref.PesticideStructureLD50 (mg/kg)Ref.
atrazine Molecules 23 02192 i0010.850
0.0689
[20]tebufenozide Molecules 23 02192 i0025000
405.40
[29]
propazine Molecules 23 02192 i0033.18
0.258
[21]imidacloprid Molecules 23 02192 i004131
10.62
[30]
simazine Molecules 23 02192 i0055
0.405
[24]acetamiprid Molecules 23 02192 i00618
14.92
[31]
carbofuran Molecules 23 02192 i0072
0.162
[25]diuron Molecules 23 02192 i008500
40.54
[22]
monocrotophos Molecules 23 02192 i00914
1.135
[26]monuron Molecules 23 02192 i0101700
437.84
[23]
dimethoate Molecules 23 02192 i01160
4.86
[27]linuron Molecules 23 02192 i0122400
194.59
[19]
carbaryl Molecules 23 02192 i013100
8.11
[28]acetylcholine Molecules 23 02192 i014
Table 2. Test set pesticides structures and their oral acute toxicities (LD50 values) against Mus musculus (upper values); test set pesticides calculated oral acute toxicities (LD50 values) against Homo sapiens (lower values, bold and italic) [32].
Table 2. Test set pesticides structures and their oral acute toxicities (LD50 values) against Mus musculus (upper values); test set pesticides calculated oral acute toxicities (LD50 values) against Homo sapiens (lower values, bold and italic) [32].
PesticideStructureLD50 (mg/kg)Ref.PesticideStructureLD50 (mg/kg)Ref.
azamethiphos Molecules 23 02192 i0151040
84.324
[44]phosmet Molecules 23 02192 i016113
9.162
[53]
azinphos-methyl Molecules 23 02192 i0177
0.567
[54]TCVP Molecules 23 02192 i018465
37.702
[55]
chlorpyrifos Molecules 23 02192 i0192000
162.162
[56]terbufos Molecules 23 02192 i0201.6
0.129
[57]
DDVP Molecules 23 02192 i02117
1.378
[58]methiocarb Molecules 23 02192 i022350
28.378
[59]
diazinon Molecules 23 02192 i02366
5.351
[60]methomyl Molecules 23 02192 i02412
0.972
[61]
fenitrothion Molecules 23 02192 i02550
40.540
[62]oxamyl Molecules 23 02192 i0265.4
0.437
[63]
glyphosate Molecules 23 02192 i0275000
405.405
[64]DDT Molecules 23 02192 i028113
9.162
[14]
malathion Molecules 23 02192 i029290
23.513
[65]2,4-D Molecules 23 02192 i030639
51.810
[66]
methyl parathion Molecules 23 02192 i03118
1.459
[67]dicamba Molecules 23 02192 i03275
61.378
[68]
naled (dibrom) Molecules 23 02192 i033160
12.972
[69]DEET Molecules 23 02192 i0341.95
0.158
[70]
parathion Molecules 23 02192 i0352
0.162
[71]sulfoxaflor Molecules 23 02192 i036750
60.810
[72]
Table 3. Training set pesticides with the known crystal structures: chemical structures, conformational analysis, superposition of the experimental conformation and generated global minima using the best performing force fields.
Table 3. Training set pesticides with the known crystal structures: chemical structures, conformational analysis, superposition of the experimental conformation and generated global minima using the best performing force fields.
PesticideFF aEglob_minbNGMS cNF dPesticideCAA f
(kJ/mol) Alignment eRMSD (Å) g
simazineMM3NA hNANA Molecules 23 02192 i038EC/MMFF
Molecules 23 02192 i037AMBER94NANANA0.667
MMFF−1022.181362782
MMFFsNANANA
OPLSAA−793.611582306
monocrotophosMM3NANANA Molecules 23 02192 i040EC/MMFF
Molecules 23 02192 i039AMBER94NANANA0.620
MMFF−307.3624162EC/MMFFs
MMFFs−310.26218613.521
OPLSAANANANA
dimethoateMM3NANANA Molecules 23 02192 i042EC/MMFF
Molecules 23 02192 i041AMBER94NANANA0.148
MMFF−383.6834044EC/MMFFs
MMFFs−386.71215403.227
OPLSAA−228.6699815EC/OPLSAA
2.746
acetamipridMM3NANANA Molecules 23 02192 i044EC/MMFF
Molecules 23 02192 i043AMBER94NANANA0.783
MMFF−31.45191614EC/MMFFs
MMFFs−31.29181660.823
OPLSAANANANA
a Force field; b Energy of the global minimum; c Number of times that a single global minimum structure was found in 10,000 processed structures; d Number of families i.e., different conformations found in 10,000 processed structures; e Red—experimental conformation, Violet—MMFF conformation, Blue—MMFFs conformation, Green—OPLSAA conformation; f Conformation dependent alignment accuracy; g RMSD measured between the heavy atoms of pesticides of experimental and best performing force field conformations; h Not available.
Table 4. Training set pesticides of the highest toxicity with the unknown crystal structures: chemical structures, conformational analysis, superposition of the experimental conformation and generated global minima using the best performing force fields.
Table 4. Training set pesticides of the highest toxicity with the unknown crystal structures: chemical structures, conformational analysis, superposition of the experimental conformation and generated global minima using the best performing force fields.
PesticideFF aEglob_min bNGMS cNF dPesticideFFDAA f
(kJ/mol) Alignment eRMSD (Å) g
atrazineMM3−1161.4844523 Molecules 23 02192 i046MM3/MMFF
Molecules 23 02192 i045AMBER94NA hNANA0.892
MMFF−1007.3214762MM3/MMFFs
MMFFs−992.63509200.931
OPLSAA−783.5736419MMFF/MMFFS
0.343
carbofuranMM329.6036172 Molecules 23 02192 i048MMFF/MMFFS
Molecules 23 02192 i047AMBER94NANANA0.075
MMFF20.3033072MMFF/OPLSAA
MMFFs15.92338420.216
OPLSAA73.27562175MMFFs/OPLSAA
0.239
a Force field; b Energy of the global minimum; c Number of times that a single global minimum structure was found in 10,000 processed structures; d Number of families i.e., different conformations found in 10,000 processed structures; e Red—experimental conformation, Violet—MMFF conformation, Blue—MMFFs conformation, Green—OPLSAA conformation; f Force field dependent alignment accuracy; g RMSD measured between the heavy atoms of pesticides experimental and best performing force field conformations; h Not available.

Share and Cite

MDPI and ACS Style

Mladenović, M.; Arsić, B.B.; Stanković, N.; Mihović, N.; Ragno, R.; Regan, A.; Milićević, J.S.; Trtić-Petrović, T.M.; Micić, R. The Targeted Pesticides as Acetylcholinesterase Inhibitors: Comprehensive Cross-Organism Molecular Modelling Studies Performed to Anticipate the Pharmacology of Harmfulness to Humans In Vitro. Molecules 2018, 23, 2192. https://doi.org/10.3390/molecules23092192

AMA Style

Mladenović M, Arsić BB, Stanković N, Mihović N, Ragno R, Regan A, Milićević JS, Trtić-Petrović TM, Micić R. The Targeted Pesticides as Acetylcholinesterase Inhibitors: Comprehensive Cross-Organism Molecular Modelling Studies Performed to Anticipate the Pharmacology of Harmfulness to Humans In Vitro. Molecules. 2018; 23(9):2192. https://doi.org/10.3390/molecules23092192

Chicago/Turabian Style

Mladenović, Milan, Biljana B. Arsić, Nevena Stanković, Nezrina Mihović, Rino Ragno, Andrew Regan, Jelena S. Milićević, Tatjana M. Trtić-Petrović, and Ružica Micić. 2018. "The Targeted Pesticides as Acetylcholinesterase Inhibitors: Comprehensive Cross-Organism Molecular Modelling Studies Performed to Anticipate the Pharmacology of Harmfulness to Humans In Vitro" Molecules 23, no. 9: 2192. https://doi.org/10.3390/molecules23092192

APA Style

Mladenović, M., Arsić, B. B., Stanković, N., Mihović, N., Ragno, R., Regan, A., Milićević, J. S., Trtić-Petrović, T. M., & Micić, R. (2018). The Targeted Pesticides as Acetylcholinesterase Inhibitors: Comprehensive Cross-Organism Molecular Modelling Studies Performed to Anticipate the Pharmacology of Harmfulness to Humans In Vitro. Molecules, 23(9), 2192. https://doi.org/10.3390/molecules23092192

Article Metrics

Back to TopTop