The Targeted Pesticides as Acetylcholinesterase Inhibitors - MDPI

26 downloads 0 Views 4MB Size Report
Aug 30, 2018 - bonding between the urea carbonyl group and mTyr124 (31% of simulation time, ...... Occupational Health Services, Inc. MSDS for Dimethoate.
molecules Article

The Targeted Pesticides as Acetylcholinesterase Inhibitors: Comprehensive Cross-Organism Molecular Modelling Studies Performed to Anticipate the Pharmacology of Harmfulness to Humans In Vitro Milan Mladenovi´c 1, *, Biljana B. Arsi´c 2,3 ID , Nevena Stankovi´c 1 , Nezrina Mihovi´c 1 , Rino Ragno 4,5 ID , Andrew Regan 6 , Jelena S. Mili´cevi´c 7 ID , Tatjana M. Trti´c-Petrovi´c 7 and Ružica Mici´c 8 1

2 3 4 5 6 7 8

*

Kragujevac Center for Computational Biochemistry, Faculty of Science, University of Kragujevac, Radoja Domanovi´ca 12, P.O. Box 60, 34000 Kragujevac, Serbia; [email protected] (N.S.); [email protected] (N.M.) Department of Mathematics, Faculty of Sciences and Mathematics, University of Niš, Višegradska 33, 18000 Niš, Serbia; [email protected] Division of Pharmacy and Optometry, University of Manchester, Oxford Road, Manchester M13 9PT, UK 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; [email protected] Alchemical Dynamics srl, 00125 Rome, Italy School of Chemistry, University of Manchester, Oxford Road, Manchester M13 9PL, UK; [email protected] Vinˇca Institute of Nuclear Sciences, University of Belgrade, P.O. Box 522, 11001 Belgrade, Serbia; [email protected] (J.S.M.); [email protected] (T.M.T.-P.) Faculty of Sciences and Mathematics, University of Priština, Lole Ribara 29, 38220 Kosovska Mitrovica, Serbia; [email protected] Correspondence: [email protected]; Tel.: +381-34-336-223

Academic Editor: Diego Muñoz-Torrero Received: 6 August 2018; Accepted: 24 August 2018; Published: 30 August 2018

 

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 LD50 s 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. Keywords: pesticides; AChE; conformational analysis; QSAR; molecular docking; molecular dynamics; quantum-chemical studies; concentration-dependent kinetic studies

Molecules 2018, 23, 2192; doi:10.3390/molecules23092192

www.mdpi.com/journal/molecules

Molecules 2018, 23, 2192

2 of 37

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–31]. Lethal doses for Homo sapiens were calculated (Table 1) according to the recommendations of Reagan-Shaw et al. [32], using Equation (1): HED = animal dose

animal Km , human Km

(1)

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].

Molecules 2018, 23, 2192

3 of 37

Molecules 2018, 23, Molecules 2018, Molecules 2018, 23, xxx23, Molecules 2018, 23, xxx Molecules 2018, 23, Molecules 2018, 23, Molecules 2018, 23, xx Molecules 2018, 23, x Molecules 2018, 23, x Molecules 2018, 23, Molecules 2018, 23, x Molecules Molecules 2018, 2018, 23, 23, xx

3 of 37 of3337 37of of 37 37 333 of 337of of 37 of of 33337 37of of 37 37 33 of 3 of 37 37

Table 1. Acetylcholine andand training set pesticides structures as well as their oraloral acute toxicities (LD(LD 50 50 Table 1. Acetylcholine training set pesticides structures as well as their acute toxicities

Table 1. Acetylcholine andand training set pesticides structures as well well as their their oraloral acute toxicities (LD(LD 50 Table 1. Acetylcholine Acetylcholine and training set pesticides structures as well well as their oral acute toxicities (LD 50 Table 1. 1. Acetylcholine and training setpesticides pesticides structures as well astheir their oral acute toxicities (LD Table Acetylcholine training set structures as as oral acute toxicities (LD 50 Table 1. training set pesticides structures as as their acute toxicities 50 50 Table 1. Acetylcholine and training set pesticides structures as well as acute toxicities (LD Table 1. Acetylcholine and training set pesticides structures as well as their oral acute toxicities (LD 50 50 Table 1.against Acetylcholine and training setvalues); pesticides structures as pesticides well as calculated their oraloral acute toxicities (LD 50 Table 1. Acetylcholine and training set pesticides structures as well as their oral acute toxicities (LD 50 values) Mus musculus (upper training set pesticides oral acute toxicities values) against Mus musculus (upper values); training set calculated oral acute toxicities values) against Mus musculus (upper values); training set pesticides calculated oral acute toxicities Table 1. Acetylcholine and training set pesticides structures as well as their oral acute toxicities (LD 50 Table 1. Acetylcholine and training set pesticides structures as well as their oral acute toxicities (LD 50 values) against Mus musculus (upper values); training set pesticides calculated oral acute toxicities Table 1. Acetylcholine and training set pesticides structures as well as their oral acute toxicities (LD 50 50 values) against Mus musculus (upper values); training set pesticides calculated oral acute toxicities values) against Mus musculus (upper values); training set pesticides calculated oral acute toxicities values) against Mus musculus (upper values); training set pesticides calculated oral acute toxicities values) against Mus musculus (upper values); training set pesticides calculated oral acute toxicities values) against Mus musculus (upper values); training set pesticides calculated oral acute toxicities values) Mus musculus (upper values); training setand pesticides calculated oraloral acute toxicities values) against Mus musculus (upper values); training set pesticides calculated oral acute toxicities (LD 50 against Homo sapiens (lower values, bold and italic) [32]. 50against values) against Homo sapiens (lower values, bold (LD(LD 50 values) values) against Homo sapiens (lower values, bold and italic) [32].[32]. values) against Mus musculus (upper values); training set pesticides calculated oral acute toxicities (LD 50 values) against Homo sapiens (lower values, bold and italic) [32]. values) against Mus musculus (upper values); training set pesticides calculated acute toxicities values) against Mus musculus (upper values); training setitalic) pesticides calculated oral acute toxicities (LD 50 values) against Homo sapiens (lower values, bold and italic) [32]. (LD 50 values) values) against Homo sapiens (lower values, bold and italic) [32]. (LD values) against Homo sapiens (lower values, bold and italic) [32]. 50 against Homo sapiens (lower values, italic) (LD 50 values) against Homo sapiens (lower values, bold and italic) [32]. 50(LD (LD 50 values) against Homo sapiens (lower values, boldbold andand italic) [32].[32]. (LD 50 values) against Homo sapiens (lower values, bold and italic) [32]. (LD(LD 50 against Homo sapiens (lower values, boldbold andand italic) [32].[32]. 50 against Homo sapiens (lower values, italic) (LD 50values) values) against Homo sapiens (lower values, bold and italic) [32]. 50 values) 50 LD LD 50 LD50 50 Ref. LD50 Pesticide Structure Pesticide Structure LD Pesticide Structure Pesticide Structure 50 LD Pesticide Structure Ref.Ref. Pesticide Structure 50 LD 50 LD Pesticide Structure (mg/kg) Ref. Pesticide Pesticide Structure Pesticide Structure Ref. Structure (mg/kg) 50 LD Pesticide Structure Ref. Pesticide Structure 50 LD (mg/kg) Pesticide Structure Ref. Pesticide Structure Pesticide Structure Pesticide Structure (mg/kg) LD 50 LD LD LD5050 Ref. 50 50 (mg/kg) Pesticide Structure Ref. Pesticide Structure (mg/kg) Pesticide Structure Ref. Pesticide Structure Pesticide Structure Ref. Pesticide Structure (mg/kg) (mg/kg) Ref. Pesticide Structure Pesticide Structure Ref. Pesticide Pesticide Structure (mg/kg) Ref. Pesticide Structure Pesticide Structure (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg)

50 LD LD 50 LD50 50 Ref. LD50 LD 50 LD Ref.Ref. 50 LD 50 LD Ref. (mg/kg) Ref. (mg/kg) 50 LD Ref. 50 LD (mg/kg) Ref. (mg/kg) LD 50 LD 50 50 Ref. LD 50 50 (mg/kg) Ref. (mg/kg) Ref. Ref. (mg/kg) (mg/kg) Ref. Ref. Ref. (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg)

atrazine atrazine atrazine atrazine atrazine atrazine atrazine atrazine atrazine atrazine atrazine atrazine atrazine atrazine

0.850 0.850 0.850 0.850 0.850 0.850 0.850 0.850 0.0689 0.0689 0.850 0.850 0.0689 0.0689 0.850 0.850 0.850 0.850 0.0689 0.0689 0.0689 0.0689 0.0689 0.0689 0.0689 0.0689 0.0689 0.0689

tebufenozide [20] tebufenozide tebufenozide [20][20] tebufenozide [20]tebufenozide [20] tebufenozide [20] tebufenozide [20] tebufenozide [20] tebufenozide [20] tebufenozide [20] [20] tebufenozide [20] tebufenozide [20]tebufenozide [20] tebufenozide

5000 5000 5000 5000 5000 5000 5000 5000 405.40 405.40 5000 5000 405.40 405.40 5000 5000 5000 405.40 405.40 405.40 405.40 405.40 405.40 405.40 405.40 405.40

[29] [29][29] [29] [29] [29] [29] [29] [29] [29] [29][29] [29] [29]

propazine propazine propazine propazine propazine propazine propazine propazine propazine propazine propazine propazine propazine propazine

3.18 3.183.18 3.18 3.18 3.18 3.18 3.18 0.258 0.258 3.18 3.18 0.258 0.258 3.18 3.18 3.18 3.18 0.258 0.258 0.258 0.258 0.258 0.258 0.258 0.258 0.258 0.258

imidacloprid [21] imidacloprid imidacloprid [21][21] imidacloprid [21]imidacloprid [21] imidacloprid [21] imidacloprid [21] imidacloprid [21] imidacloprid [21] imidacloprid [21] [21] imidacloprid [21] imidacloprid [21]imidacloprid [21] imidacloprid

131 131131 131 131 131 131 131 10.62 13110.62 131 10.62 10.62 131 131 131 10.62 10.62 10.62 10.62 10.62 10.62 10.62 10.62 10.62

[30] [30][30] [30] [30] [30] [30] [30] [30] [30] [30][30] [30] [30]

simazine simazine simazine simazine simazine simazine simazine simazine simazine simazine simazine simazine simazine

5 5 555 555 0.405 0.405 0.405 5550.405 555 0.405 0.405 0.405 0.405 0.405 0.405 0.405 0.405 0.405 0.405

acetamiprid [24] acetamiprid acetamiprid [24][24] acetamiprid [24]acetamiprid [24] acetamiprid [24] acetamiprid [24] acetamiprid [24] acetamiprid [24] acetamiprid [24] [24] acetamiprid [24] acetamiprid [24] acetamiprid [24]acetamiprid

18 18 18 18 18 18 18 18 14.92 18 18 14.92 14.92 1814.92 18 18 14.92 14.92 14.92 14.92 14.92 14.92 14.92 14.92 14.92

[31] [31][31] [31] [31] [31] [31] [31][31] [31] [31] [31] [31] [31]

2 2 2222 222 0.162 0.162 0.162 22 0.162 222 0.162 0.162 0.162 0.162 0.162 0.162 0.162 0.162 0.162 0.162

[25] [25][25] [25] [25] [25] [25] [25] [25] [25] [25] [25][25] [25]

diuron diuron diuron diuron diuron diuron diuron diuron diuron diuron diuron diuron diuron diuron

500 500500 500 500 500 500 500 40.54 500 500 40.54 40.54 50040.54 500 500 40.54 40.54 40.54 40.54 40.54 40.54 40.54 40.54 40.54

[22] [22][22] [22] [22] [22] [22] [22] [22] [22] [22] [22][22] [22]

monocrotophos monocrotophos monocrotophos monocrotophos monocrotophos monocrotophos monocrotophos monocrotophos monocrotophos monocrotophos monocrotophos monocrotophos monocrotophos monocrotophos

14 14 14 14 14 14 14 14 14 1.135 141.135 14 1.135 1.135 14 14 14 1.135 1.135 1.135 1.135 1.135 1.135 1.135 1.135 1.135 1.135

[26] [26] [26] [26] [26] [26] [26] [26] [26] [26] [26] [26][26] [26]

monuron monuron monuron monuron monuron monuron monuron monuron monuron monuron monuron monuron monuron monuron

1700 1700 1700 1700 1700 1700 1700 1700 437.84 437.84 1700 1700 437.84 437.84 1700 1700 1700 437.84 437.84 437.84 437.84 437.84 437.84 437.84 437.84 437.84

[23] [23][23] [23] [23] [23] [23] [23] [23] [23] [23] [23][23] [23]

dimethoate dimethoate dimethoate dimethoate dimethoate dimethoate dimethoate dimethoate dimethoate dimethoate dimethoate dimethoate dimethoate

60 60 60 60 60 60 60 60 60 4.86 60 4.86 60 4.86 4.86 4.86 60 60 60 4.86 4.86 4.86 4.864.86 4.86 4.86 4.86 4.86

[27] [27] [27] [27] [27] [27] [27] [27] [27][27] [27] [27] [27] [27]

linuron linuron linuron linuron linuron linuron linuron linuron linuron linuron linuron linuron linuron linuron

2400 2400 2400 2400 2400 2400 2400 2400 194.59 2400 194.59 2400 194.59 194.59 2400 2400 2400 194.59 194.59 194.59 194.59 194.59 194.59 194.59 194.59 194.59

[19] [19][19] [19] [19] [19] [19] [19] [19][19] [19] [19] [19] [19]

carbaryl carbaryl carbaryl carbaryl carbaryl carbaryl carbaryl carbaryl carbaryl carbaryl carbaryl carbaryl carbaryl

100 100100 100 100 100 100 100 100 8.11 8.11 8.11 100 100 8.11 8.11 100 100 100 8.11 8.11 8.11 8.118.11 8.11 8.11 8.11 8.11

acetylcholine [28] acetylcholine [28] acetylcholine [28] acetylcholine [28] acetylcholine [28]acetylcholine [28] acetylcholine [28] acetylcholine acetylcholine [28] acetylcholine [28][28] acetylcholine [28] acetylcholine [28] [28] acetylcholine [28] acetylcholine

carbofuran carbofuran carbofuran carbofuran carbofuran carbofuran carbofuran carbofuran carbofuran carbofuran carbofuran carbofuran carbofuran

O O O O O O O O O NO O O N N OO O HN O OO H N O N H N H N O O O H H N H O N O O H N O N O O H N H O H HH O O

O O O O O O O

For either mouse or humans as aas model organism, AChE is homodimer (Figure 1a), formed of For either mouse or model organism, AChE aa homodimer (Figure 1a), formed of For either mouse or humans asas model organism, AChE is aaais homodimer (Figure 1a), formed offormed For either mouse or humans humans as model organism, AChE is homodimer (Figure 1a), formed of For either mouse orhumans humans aaaaamodel organism, AChE a homodimer (Figure 1a), For either mouse or as aaaas model organism, AChE is homodimer (Figure 1a), formed of For either mouse or humans as model organism, AChE is homodimer (Figure 1a), formed of For either mouse or humans model organism, AChE is aaaishomodimer (Figure 1a), formed of For either mouse or humans as model organism, AChE is a homodimer (Figure 1a), formed of For either mouse or humans as a model organism, AChE is a homodimer (Figure 1a), formed of For either mouse or humans as a model organism, AChE is homodimer (Figure 1a), formed of two independent units, comprising 614 amino acids [36]. The AChE active site comprises the aromatic two independent units, comprising 614 amino acids [36]. The AChE active site comprises the aromatic two independent units, comprising 614 amino acids [36]. The AChE active site comprises the aromatic For either mouse or humans as a model organism, AChE is a homodimer (Figure 1a), formed of For either mouse or humans as a model organism, AChE is a homodimer (Figure 1a), formed of For either mouse or humans as a model organism, AChE is a homodimer (Figure 1a), formed of the two independent units, comprising 614 amino acids [36]. The AChE active site comprises the aromatic independent units, comprising 614 amino acids [36]. The AChE active site comprises the aromatic two independent units, comprising 614 amino acids [36]. The AChE active site comprises aromatic of two two independent units, comprising 614 amino acids [36]. The AChE active sitethe comprises two independent units, comprising amino acids [36]. The AChE active site comprises aromatic two independent units, comprising 614 amino acids [36]. The AChE active site comprises the aromatic two independent units, comprising 614614 amino acids [36]. The AChE active site comprises thethe aromatic 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 anionic domain, which accommodates the ACh positive quaternary amine, and the esterase catalytic anionic domain, which accommodates the ACh positive quaternary amine, and the esterase catalytic two independent units, comprising 614 amino acids [36]. The AChE active site comprises the aromatic two independent units, comprising 614 amino acids [36]. The AChE active site comprises the aromatic 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 anionic domain, which accommodates the ACh positive quaternary amine, and the esterase catalytic anionic domain, which accommodates the ACh positive quaternary amine, and the esterase catalytic aromatic anionic domain, which accommodates the ACh positive quaternary and the esterase anionic domain, which accommodates the ACh positive quaternary amine, and the esterase catalytic anionic domain, which accommodates the ACh positive quaternary amine, and the esterase catalytic anionic domain, which accommodates theSer203, ACh positive quaternary amine, and theamine, esterase catalytic 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 domain, containing the catalytic triad His447, and Glu334 with the purpose of domain, containing theaccommodates catalytic triad Ser203, His447, and Glu334 with theand purpose ofesterase hydrolysing anionic domain, which the ACh positive quaternary amine, the esterase catalytic anionic domain, which accommodates the ACh positive quaternary amine, and the catalytic 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 hydrolysing domain, containing the catalytic triad Ser203, His447, and Glu334 with the purpose of hydrolysing domain, containing the catalytic triad Ser203, His447, and Glu334 with the purpose of hydrolysing domain, containing the catalytic triad Ser203, His447, and Glu334 with the purpose of hydrolysing domain, containing the catalytic triad Ser203, His447, and Glu334 with the purpose of hydrolysing catalytic domain, containing the catalytic triad Ser203, His447, and Glu334 with the purpose of domain, containing the catalytic triad Ser203, His447, and Glu334 with the purpose ofan hydrolysing 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 acetylatedACh to acetate and choline. The ester hydrolysis reaction leads to formation of an acetylatedthethe ACh to acetate and choline. The ester hydrolysis reaction leads to the formation ofof an acetylateddomain, containing the catalytic triad Ser203, His447, and Glu334 with the purpose hydrolysing domain, containing the catalytic triad Ser203, His447, and Glu334 with the purpose of 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 anhydrolysing acetylatedthe ACh to acetate and choline. The ester hydrolysis reaction leads to the formation of an acetylatedthe ACh to acetate and choline. The ester hydrolysis reaction leads to the formation of an acetylatedthe ACh to acetate and choline. The ester hydrolysis reaction leads to the formation of an acetylatedthe ACh to acetate and choline. The ester hydrolysis reaction leads to the formation of an acetylatedhydrolysing the ACh to acetate and choline. The ester hydrolysis reaction leads to the formation of the ACh to acetate and choline. The ester hydrolysis reaction leads to the formation of an acetylatedthe ACh to acetate and choline. The ester hydrolysis reaction leads to the formation of an acetylatedenzyme. Then, the acetyl-enzyme undergoes nucleophilic attack by aby water molecule, assisted by the enzyme. Then, acetyl-enzyme undergoes nucleophilic attack aa water molecule, by enzyme. Then, thethe acetyl-enzyme undergoes nucleophilic attack by aby water molecule, assisted by the the the ACh to acetate and choline. The ester hydrolysis reaction leads to the formation ofassisted an acetylatedthe ACh to and choline. The ester hydrolysis reaction leads to the formation of an the ACh to acetate acetate and choline. The ester hydrolysis reaction leads to the formation ofassisted an acetylatedacetylatedenzyme. Then, the acetyl-enzyme undergoes nucleophilic attack by water molecule, assisted by the the enzyme. Then, the acetyl-enzyme undergoes nucleophilic attack by aaby water molecule, by enzyme. Then, the acetyl-enzyme undergoes nucleophilic attack a water water molecule, assisted by the enzyme. Then, the acetyl-enzyme undergoes nucleophilic attack a molecule, assisted by the enzyme. Then, the acetyl-enzyme undergoes nucleophilic attack by water molecule, assisted by the enzyme. Then, the acetyl-enzyme undergoes nucleophilic attack by a water molecule, assisted by the enzyme. Then, the acetyl-enzyme undergoes nucleophilic attack by a water molecule, assisted by the anenzyme. acetylated-enzyme. Then, the acetyl-enzyme undergoes nucleophilic attack by a water molecule, His447, freeing the acetic acid, likely as an acetate ion, and regenerating the free enzyme (Figure 1b) His447, freeing the acetic acid, likely as an acetate and regenerating the free enzyme (Figure His447, freeing thethe acetic acid, likely as an acetate ion,ion, andattack regenerating themolecule, free enzyme (Figure 1b) Then, the acetyl-enzyme undergoes by aby water assisted by the enzyme. Then, acetyl-enzyme undergoes nucleophilic attack aa water molecule, assisted by the enzyme. Then, the acetyl-enzyme undergoes nucleophilic attack by water molecule, assisted by 1b) the His447, freeing the acetic acid, likely as acetate annucleophilic acetate ion, and regenerating the free enzyme (Figure 1b) His447, freeing the acetic acid, likely as an ion, and regenerating the free enzyme (Figure 1b) His447, freeing the acetic acid, likely as an acetate ion, and regenerating the free enzyme (Figure 1b) His447, freeing acetic acid, likely as an acetate ion, and regenerating free enzyme (Figure 1b) His447, freeing the acetic acid, likely as an acetate ion, and regenerating the free enzyme (Figure 1b) His447, freeing thethe acetic acid, likely as an an acetate ion,as and regenerating thethe free enzyme (Figure 1b)enzyme His447, freeing the acetic acid, likely as an acetate ion, and regenerating the free enzyme (Figure 1b) [33]. [33]. assisted by the His447, freeing the acetic acid, likely an acetate ion, and regenerating the free [33]. His447, freeing the acetic acid, likely as acetate ion, and regenerating the free enzyme (Figure 1b) His447, freeing the acetic acid, likely as an acetate ion, and regenerating the free enzyme (Figure 1b) His447, freeing the acetic acid, likely as an acetate ion, and regenerating the free enzyme (Figure 1b) [33]. [33]. [33]. [33]. [33]. [33]. [33]. Given that the pharmacology of the majority of widely used pesticides (Table 1) is still poorly Given that the pharmacology of the majority of widely used pesticides (Table 1) is still poorly Given that the pharmacology of the majority of widely used pesticides (Table 1) is still poorly (Figure 1b) [33]. [33]. [33]. [33]. Given that the pharmacology of the majority of widely used pesticides (Table 1) is still poorly Given that the pharmacology of the majority of widely used pesticides (Table 1) is poorly Given that the pharmacology of the the majority of widely widely used pesticides (Table 1)still is still still poorly Given that the pharmacology of majority of used pesticides (Table is poorly Given that the pharmacology of the majority of widely used pesticides (Table 1) is still poorly Given that the pharmacology ofwas the majority ofthe widely used pesticides (Table 1) is is1) still poorly Given that the pharmacology of the majority of widely used pesticides (Table 1) is still poorly understood, the focus of this study to provide LB or SB cheminformatics-based guidelines understood, the focus of study was to provide LB or SB cheminformatics-based guidelines understood, the focus of this study was to provide the LB orused SB cheminformatics-based guidelines Given that the pharmacology of the majority of widely pesticides (Table 1) poorly Given that the pharmacology of the majority of widely pesticides (Table 1) is poorly Given that the pharmacology of the majority ofthe widely used pesticides (Table 1)still is1)still still poorly understood, focus of this this study was to provide the LB orused SB cheminformatics-based guidelines Given that the pharmacology of the majority widely used pesticides (Table is still poorly understood, the focus of this study was to provide the LB or SB cheminformatics-based guidelines understood, the focus of this study was to provide the LB or SB cheminformatics-based guidelines understood, the focus of this study was to provide the LB or SB cheminformatics-based guidelines understood, the focus of this study was to provide the LB or SB cheminformatics-based guidelines understood, the focus of this study was to provide the LB or SB cheminformatics-based guidelines 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 for several important tasks, and unresolved questions: (1) how to define the pesticides’ structures for several important tasks, and unresolved questions: (1) how to define the pesticides’ structures in in understood, the focus of this study was to provide the LB or SB cheminformatics-based guidelines understood, the focus of this study was to provide the LB or SB cheminformatics-based guidelines 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 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 for several important tasks, and unresolved questions: (1) how to define the pesticides’ structures in for several important tasks, and unresolved questions: (1) how to define the pesticides’ structures in for several important tasks, and unresolved questions: (1) how to define the pesticides’ structures in for several important tasks, and unresolved questions: (1) how to define the pesticides’ structures in 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 aqueous solution prior to interaction with AChE? (2) to define the pesticides’ bioactive aqueous solution prior to interaction with AChE? (2) how to define the pesticides’ bioactive for several important tasks, and unresolved questions: (1) how to define the pesticides’ structures in for several important tasks, and unresolved questions: (1) how to define the structures in for several important tasks, and unresolved questions: (1) how to define the pesticides’ structures in aqueous solution prior to interaction with AChE? (2) toto define the pesticides’ bioactive solution prior to interaction with AChE? (2) how to define the pesticides’ bioactive aqueous solution prior to interaction with AChE? (2) how to define the pesticides’ bioactive foraqueous several important tasks, and unresolved questions: (1) how define the pesticides’ structures aqueous solution prior to interaction with AChE? (2) how to define the pesticides’ bioactive aqueous solution prior to interaction with AChE? (2) how to define the pesticides’ bioactive aqueous solution prior to interaction with AChE? (2) how to define the pesticides’ bioactive 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 conformations upon the interaction with AChE? Other important assignments of the LB or SB conformations upon the interaction with AChE? Other important assignments of the LB or SB aqueous solution prior to interaction with AChE? (2) how to define the pesticides’ bioactive aqueous solution prior to interaction with AChE? (2) how to define the pesticides’ bioactive aqueous solution prior to interaction with AChE? (2) how to assignments define the pesticides’ bioactive conformations upon the interaction with AChE? Other important assignments of the the LB or SB SB upon the interaction with AChE? Other important the LB or SB conformations upon the interaction with AChE? Other important assignments of LB or in conformations aqueous solution prior to interaction with AChE? (2)important how to define theof pesticides’ bioactive conformations upon the interaction with AChE? Other assignments of the LB or SB conformations upon the interaction with AChE? Other important assignments of the LB or SB conformations upon the interaction with AChE? Other important assignments of(1) the LB orany SB 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: is there protocols applied herein were to answers to important questions: is protocols applied herein were to provide provide answers to several several important questions: (1)the isthe there any conformations upon the interaction with AChE? Other important assignments of LB or SB conformations upon the interaction with AChE? Other important assignments of LB or SB conformations upon the interaction with AChE? Other important assignments of(1) the LB orany SB protocols applied herein were to provide provide answers to several several important questions: (1) is there there any protocols applied herein were to answers to important questions: (1) is there any protocols applied herein were to provide answers to several important questions: (1) is there any protocols applied herein were to provide answers to several important questions: (1) is any protocols applied herein were to provide answers to several important questions: (1) is there any conformations upon the interaction with AChE? Other important assignments ofbefore the orthere SB protocols protocols applied herein were to provide provide answers to from several important questions: (1) isLB there any 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, the actual possibility to the pesticides acute toxicity pre-bound conformations, before actual possibility to applied anticipate the pesticides acute toxicity from pre-bound conformations, before the actual protocols applied herein were to answers to several important questions: (1) is there any protocols herein were to answers to several important questions: (1) is there any protocols applied herein were to provide provide answers topre-bound several important questions: (1) is the there any possibility to anticipate anticipate the pesticides acute toxicity from pre-bound conformations, before the actual possibility to anticipate the pesticides acute toxicity from conformations, before the actual possibility to anticipate the pesticides acute toxicity from pre-bound conformations, before the actual possibility to anticipate the pesticides acute toxicity from pre-bound conformations, before the actual possibility to anticipate the pesticides acute toxicity from pre-bound conformations, before the actual applied herein were toAChE? provide answers topesticides’ several important questions: (1) is from there any possibility to possibility towith anticipate the(2) pesticides acute toxicity from pre-bound conformations, before the actual possibility to anticipate the pesticides acute toxicity from pre-bound conformations, before the actual interaction AChE? can the pesticides’ acute toxicity be quantified the bound interaction with (2) can the acute toxicity be quantified from the bound interaction with AChE? (2) can theacute pesticides’ acute toxicity be quantified from the bound possibility to anticipate the pesticides toxicity from pre-bound conformations, before the actual possibility to anticipate the pesticides acute toxicity from pre-bound conformations, before the actual 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 interaction with AChE? (2) can the pesticides’ acute toxicity be quantified from the bound interaction with AChE? (2) can the pesticides’ acute toxicity be quantified from the bound interaction with AChE? (2) can the pesticides’ acute toxicity be quantified from the bound interaction with AChE? (2) can the pesticides’ acute toxicity be quantified from the bound interaction with AChE? (2) (2) can thethe pesticides’ acute toxicity be be quantified from thethe bound interaction with AChE? (2) can the pesticides’ acute toxicity be quantified from the bound anticipate thewith pesticides acute toxicity from pre-bound conformations, before the actual interaction interaction AChE? (2) can the pesticides’ acute toxicity be quantified from the bound interaction with AChE? can pesticides’ acute toxicity quantified from bound interaction with AChE? (2) can the pesticides’ acute toxicity be quantified from the bound

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

Molecules2018, 2018,23, 23,2192 x Molecules

of 37 37 44 of

conformations, upon the actual interaction with AChE? To the best of our knowledge, there are no either LB or SB efforts to predict of investigated pesticides 1) for previous reports on either LB orthe SBbioactive efforts toconformations predict the bioactive conformations of (Table investigated both mouse and human AChEs, as well as to enumerate their acute toxicities by means of free energy pesticides (Table 1) for both mouse and human AChEs, as well as to enumerate their acute toxicities of (Eglob_min ) in aqueous free energy of and/or bindingfree (∆G tobinding AChE. byformation means of free energy ofthe formation (Esolution glob_min) inand/or the aqueous solution energy binding )of All cheminformatics findings were externally validated on the series of other commonly (∆Gthe binding ) to AChE. All the cheminformatics findings were externally validated on the series of used other pesticides (Table 2), also known as AChE inhibitors (i.e., the test set). commonly used pesticides (Table 2), also known as AChE inhibitors (i.e., the test set).

Figure structure of of hAChE hAChE (PDB (PDBID: ID:4EY5) 4EY5)in incomplex complexwith withhuperzine huperzineAA(a). (a).Subunit Subunit Figure 1. 1. The The crystal crystal structure A A is depicted in pink ribbons whereas the B subunit is coloured in orange. Huperzine A is presented is depicted in pink ribbons whereas the B subunit is coloured in orange. Huperzine A is presented in in black and structureisishighlighted highlightedwith withdashed dashedlines. lines. The The anionic anionic sub-site sub-site is is described black and itsitsstructure described by by blue blue spheres, while the esterase sub-area is encircled by red spheres. For the sake of clarity, hydrogen spheres, while the esterase sub-area is encircled by red spheres. For the sake of clarity, hydrogen atoms atoms were were omitted. omitted. The The graphic graphic was was generated generated using using the the UCSF UCSF Chimera Chimera software software (Resource (Resource for for Biocomputing, Visualization, and Informatics (RBVI), University of California, San Francisco, CA, Biocomputing, Visualization, and Informatics (RBVI), University of California, San Francisco, CA, USA). USA). Biochemical Biochemical mechanism mechanismof ofACh AChhydrolysis hydrolysisby byAChE AChE(b). (b).

The test test set set was was compiled compiled of of organophosphates, organophosphates, carbamates carbamates and and structurally structurally non-classified non-classified The pesticides. Regarding organophosphates, chlorpyrifos [37], DDVP [14], fenitrothion [38], azinphospesticides. methyl [39], [39], naled naled (dibrom) (dibrom) [40], [40], TCVP TCVP [41], [41], parathion [39], methyl parathion [42], diazinon [43], methyl phosmet [14], [14], azamethiphos azamethiphos [44], [44], and and terbufos terbufos [45], [45], respectively, respectively, are are used used as as insecticides. insecticides. Moreover, phosmet glyphosate is is aa broad-spectrum broad-spectrum systemic systemic herbicide herbicide and and crop crop desiccant desiccant [46], [46], malathion malathion is is mostly mostly used used glyphosate for mosquito mosquito eradication eradication [47], [47], parathion parathion isis used used as as an anacaricide. acaricide. Phosmet Phosmet is is mainly mainly used used on on apple apple for trees to to control control codling codling moth moth and and on on aa numerous numerous fruit fruit crops, crops, ornamentals, ornamentals, and and wines wines for for the the control control trees of aphids, aphids, suckers, suckers, mites, mites, and and fruit fruit flies. flies. Terbufos Terbufosisisaa nematicide nematicideused used on on corn, corn, sugar sugar beets beets and and grain grain of sorghum,and andit it is extremely to birds. thehand, othercarbamate hand, carbamate such as sorghum, is extremely toxictoxic to birds. On theOn other pesticides,pesticides, such as methomyl methomyl and methiocarb are used against insects; oxamyl is an extremely toxicto pesticide tofish, humans, and methiocarb are used against insects; oxamyl is an extremely toxic pesticide humans, and fish, and birds [48–50], whereas methiocarb is additionally used as a bird repellent, acaricide and

Molecules 2018, 23, 2192

5 of 37

birds [48–50], whereas methiocarb is additionally used as a bird repellent, acaricide and molluscicide. Molecules 2018, x 37 Molecules 2018, 23, 23, x 23, 5 of of 37 37 555 of Molecules 2018, 23, of 37 Molecules 2018, Molecules 2018, 23, xxxx of 37 37 Molecules 2018, Molecules 2018, 23, of 37 Molecules 2018, 23, xxx 5555 of 37 Molecules 2018, 23, of 37 Molecules 2018, 23, x 23, ofmimics 37 55 of Among the structurally non-classified pesticides, 2,4-D is widely used as a herbicide, which the action of the plant growth hormone auxin, resulting in the uncontrolled growth and eventual death in molluscicide. Among the structurally non-classified pesticides, 2,4-D is widely used as herbicide, molluscicide. Among the structurally non-classified pesticides, 2,4-D is widely used as herbicide, molluscicide. Among the structurally non-classified pesticides, 2,4-D is widely used as aaa herbicide, molluscicide. Among the non-classified pesticides, 2,4-D is used as molluscicide. Among the structurally structurally non-classified pesticides, 2,4-D is widely widely used as aaaaa herbicide, herbicide, molluscicide. Among the non-classified pesticides, is as molluscicide. Among the structurally non-classified pesticides, 2,4-D is widely used as herbicide, molluscicide. Among the structurally non-classified pesticides, 2,4-D is widely used as herbicide, molluscicide. Among the structurally structurally non-classified pesticides, 2,4-D2,4-D is widely widely used used as a herbicide, herbicide, susceptible plants [14]; DDT is originally developed as an insecticide, used to control malaria and typhus, which mimics the action of the plant growth hormone auxin, resulting in the uncontrolled growth which mimics the action of the plant growth hormone auxin, resulting in the uncontrolled growth which mimics the action of the plant growth hormone auxin, resulting in the uncontrolled growth which mimics the of growth hormone auxin, resulting in uncontrolled growth which mimics the action action ofplant the plant plant growth hormone auxin, resulting in the uncontrolled uncontrolled growth mimics the the growth hormone resulting the growth which mimics the action of the plant growth hormone auxin, resulting in the uncontrolled growth which mimics the action of the plant growth hormone auxin, resulting in the uncontrolled growth whichwhich mimics the action action of the theof plant growth hormone auxin,auxin, resulting in the thein uncontrolled growth and eventual death in susceptible plants [14]; DDT is originally developed as an insecticide, used to and eventual death in susceptible plants [14]; DDT is originally developed as an insecticide, used to and as a contact poison against several arthropods [51]; DEET is the most common active ingredient and eventual death in susceptible plants [14]; DDT is originally developed as an insecticide, used to and eventual death in susceptible plants [14]; DDT is originally developed as an insecticide, used to and eventual eventual death in susceptible susceptible plants [14]; is DDT is originally originally developed as an insecticide, insecticide, used to to and in [14]; DDT is developed an used and death in plants [14]; originally developed as insecticide, used and eventual death in susceptible plants [14]; DDT is originally developed as an insecticide, used to and eventual eventual deathdeath in susceptible susceptible plantsplants [14]; DDT DDT is originally developed as an an as insecticide, used to to control malaria and typhus, and as a contact poison against several arthropods [51]; DEET is the most control malaria and typhus, and as a contact poison against several arthropods [51]; DEET is the most control malaria and typhus, and as a contact poison against several arthropods [51]; DEET is the most incontrol insect repellents for protection against mosquitoes, ticks, fleas, chiggers, leeches and many biting control malaria and typhus, and as a contact poison against several arthropods [51]; DEET is the most control malaria and typhus, and as a contact poison against several arthropods [51]; DEET is the most control malaria and typhus, and as a contact poison against several arthropods [51]; DEET is the most malaria and and aa contact poison against several arthropods [51]; is control malaria and typhus, as a contact poison against several arthropods [51]; DEET is the most control malaria and typhus, typhus, and as asand contact poison against several arthropods [51]; DEET DEET is the the most most common active ingredient in insect repellents for protection against mosquitoes, ticks, fleas, chiggers, common activeactive ingredient in insect insect repellents for protection protection against mosquitoes, ticks, ticks, fleas, chiggers, common active ingredient in insect repellents for protection against mosquitoes, ticks, fleas, chiggers, common active ingredient in repellents for against mosquitoes, ticks, fleas, chiggers, common active ingredient in insect insect repellents for protection protection against mosquitoes, ticks, fleas, chiggers, insects [52]. common ingredient in repellents for against mosquitoes, chiggers, common active ingredient in repellents for against mosquitoes, ticks, chiggers, common ingredient in insect repellents for protection against mosquitoes, ticks, fleas, chiggers, common activeactive ingredient in insect insect repellents for protection protection against mosquitoes, ticks, fleas, fleas, fleas, chiggers, leeches and many biting insects [52]. leeches and many many bitingbiting insects [52]. [52]. leeches and many biting insects [52]. leeches and biting insects [52]. leeches and many many biting insects [52]. leeches and insects leeches and biting insects [52]. leeches and many insects leeches and many many bitingbiting insects [52]. [52]. Table 2. Test set pesticides structures and their oral acute toxicities (LD50 values) against Mus musculus Table 2. Test set pesticides structures and their oral acute toxicities (LD 50 against Mus musculus 50 values) 50 values) Table 2. Test set pesticides structures and their oral acute toxicities (LD 50 against Mus musculus Table 2. Test set pesticides structures and their oral acute toxicities (LD 50 values) against Mus musculus 50 50 Table 2. Test set pesticides structures and their oral acute toxicities (LD values) against Mus musculus Table 2. Test setset pesticides structures andoral their oral acute toxicities (LD50 values) against Mus musculus 50(LD 50 values) Test set pesticides structures and their oral acute toxicities (LD 50 against Mus musculus Table 2. Test set pesticides structures and their oral acute toxicities (LD 50 values) against Mus musculus (upper values); test pesticides calculated oral acute toxicities values) against Homo sapiens Table 2. Test set pesticides structures and acute toxicities (LD against Mus musculus Table 2. Test set pesticides structures and their oral acute toxicities (LD 50 values) against Mus musculus TableTable 2. Test2. set pesticides structures and their their oral acute toxicities (LD50 50 values) values) against Mus musculus (upper values); test set pesticides calculated oral acute toxicities (LD 50 values) against Homo sapiens 50 50 (upper values); test set pesticides calculated oral acute toxicities (LD 50 values) against Homo sapiens (upper values); test set pesticides calculated oral acute toxicities (LD 50 values) against Homo sapiens 50 50 (upper values); test set pesticides calculated oral acute toxicities (LD values) against Homo sapiens (upper values); test set pesticides calculated oral acute toxicities (LD values) against Homo sapiens 50 50 (upper values); test set pesticides calculated oral acute toxicities (LD 50 values) against Homo sapiens (upper values); test set pesticides calculated oral acute toxicities (LD 50 values) against Homo sapiens (upper values); test set pesticides calculated oral acute toxicities (LD 50 values) against Homo sapiens (upper values); test set pesticides calculated oral acute toxicities (LD 50 values) against Homo sapiens (lower andpesticides italic) [32]. (uppervalues, values);bold test set calculated oral acute toxicities (LD50 values) against Homo sapiens

(lower values, bold and italic) [32]. (lower(lower values,values, bold and and italic) [32]. (lower values, bold and italic) [32]. (lower values, bold italic) [32]. (lower values, bold and italic) italic) [32]. [32]. bold and (lower values, bold and italic) [32]. (lower values, bold italic) [32]. bold and italic) [32]. (lower(lower values,values, bold and and italic) [32]. 50 50 50 LD 50 LD 50 50 LD50 LD50 50 50 50 50 LD LD LD50 50 50 50 LD LD 50 LD LD 5050 LD 50 LD 50Ref. 50Ref. Ref. Pesticide Structure Ref. Pesticide Structure LD LD Pesticide Structure Pesticide Structure Ref.Ref. Pesticide Ref. Pesticide Pesticide Structure Ref. Pesticide Structure Ref. 50 LD50Ref. 50 LD50Ref. LD50 LD50 Pesticide Structure Pesticide Structure Pesticide Structure Structure (mg/kg) Ref. Pesticide Structure Structure (mg/kg) Ref. Pesticide Structure Pesticide Structure (mg/kg) (mg/kg) Pesticide Structure Ref. Pesticide Structure Ref. Pesticide Ref. Pesticide Ref. (mg/kg) (mg/kg) Pesticide Structure Ref. Pesticide Structure Ref. (mg/kg) (mg/kg) Pesticide Structure Ref. Pesticide Structure Ref. Ref. (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) (mg/kg) 113 1040 113 113 1040 113 1040 113 1040 1040 1040 113 [53] [53] 1040 [44] 113 azamethiphos phosmet [44] 113 1040 113 1040 azamethiphos phosmet azamethiphos phosmet [53] [44] 113 1040 113 9.162 1040 84.324 azamethiphos [44] [44] phosmet [53] phosmet [53] [44] azamethiphos phosmet [53] [44] azamethiphos phosmet [53] [44] azamethiphos phosmet azamethiphos azamethiphos phosmet [53] [44] 9.162 84.324 azamethiphos phosmet [53] [44] 84.324 9.162 9.162 84.324 azamethiphos phosmet [53] [53] [44] 9.162 84.324 9.162 84.324 9.162 84.324 9.162 84.324 9.162 84.324 9.162 9.162 84.32484.324

azinphosazinphosazinphosazinphosazinphosazinphosazinphosazinphos-methyl azinphosazinphosmethyl methyl methyl methyl methyl methyl methyl methyl methyl methyl chlorpyrifos chlorpyrifos chlorpyrifos chlorpyrifos chlorpyrifos chlorpyrifos chlorpyrifos chlorpyrifos chlorpyrifos chlorpyrifos chlorpyrifos

465 465 465 465 465 465 465 [55] 465 [55] [55] [55] 465 [55] 465 37.702 [55] [55] [55] 37.702 [55] 37.702 37.702 [55] [55] 37.702 37.702 37.702 37.702 37.702 37.702 37.702

7 7777 7777 [54] TCVP TCVP [54] [54] [54] TCVP TCVP [54] 7 0.567 TCVP TCVP [54]TCVP TCVP [54] 0.567 [54] TCVP [54] 0.567 0.567 TCVPTCVP [54] [54] 0.567 0.567 0.567 0.567 0.567 0.567 0.567 2000 2000 2000 2000 2000 2000 2000 [56] 2000 terbufos [56] [56] terbufos 2000 terbufos 2000 terbufos [56] 2000162.162 terbufos [56] terbufos [56] terbufos 162.162 terbufos [56] terbufos [56] terbufos [56] 162.162 162.162 terbufos [56] [56] 162.162 162.162 162.162 162.162 162.162 162.162 162.162

1.6 1.6 1.6 1.6 1.6 0.129 0.129 0.129 0.129 0.129

1.6 1.6 1.6 1.6 1.6 [57] [57] [57] [57] 1.6 [57] [57] [57] 0.129[57] 0.129 0.129 [57] [57] 0.129[57] 0.129 0.129

DDVP DDVPDDVP DDVP DDVP DDVP DDVP DDVP DDVP DDVPDDVP

17 1717 17 17 17 17 [58] methiocarb 17 17 methiocarb 17 [58] [58] methiocarb methiocarb [58] 17 1.378 methiocarb [58] methiocarb [58] methiocarb methiocarb [58] methiocarb [58] 1.378 1.378 methiocarb 1.378 [58] methiocarb [58] [58] 1.378 1.378 1.378 1.378 1.378 1.378 1.378

350 350 350 350 350 350 350 [59] 350 350 [59] [59] 350 [59] [59] 350 28.378 [59] [59] [59] [59] 28.378 28.378 [59] [59] 28.378 28.378 28.378 28.378 28.378 28.37828.378

diazinon diazinon diazinon diazinon diazinon diazinon diazinon diazinon diazinon diazinon diazinon

66 66 66 66 66 [60] [60] 66 66 methomyl 66 methomyl 66 66 methomyl 66 methomyl [60] [60] [60] methomyl methomyl [60] methomyl [60] methomyl 5.351 methomyl [60] 5.351 methomyl [60] 5.351 methomyl [60] [60] 5.351 5.351 5.351 5.351 5.351 5.351 5.351 5.351

12 12 12 12 12 [61] [61] 12 12 12 12 [61] 12 0.972 [61] [61] [61] [61] [61] 0.972 0.972 [61] 0.972 [61] [61] 0.972 0.972 0.972 0.972 0.972 0.972 0.972

fenitrothion fenitrothion fenitrothion fenitrothion fenitrothion fenitrothion fenitrothion fenitrothion fenitrothion fenitrothion fenitrothion

50 50 50 50 50 [62] [62] 50 50 oxamyl 50 oxamyl 50 50 oxamyl [62] 50 oxamyl [62] oxamyl [62] [62] [62] oxamyl oxamyl [62] oxamyl 40.540 oxamyl [62] 40.540 oxamyl [62] 40.540 oxamyl [62] 40.540 40.540 40.540 40.540 40.540 40.540 40.540 40.540

5.4 5.4 5.4 5.4 5.4 0.437 0.437 0.437 0.437 0.437

5.4 5.4 5.4 [63] [63] 5.4 5.4 [63] [63] [63] [63] [63] 0.437 [63] 0.437 [63] [63] 0.437[63] 0.437 0.437

glyphosate glyphosate glyphosate glyphosate glyphosate glyphosate glyphosate glyphosate glyphosate glyphosate glyphosate

5000 5000 5000 5000 5000[64] [64] 5000 5000 DDT 5000 5000 5000 DDT DDT [64] 5000 DDT [64] DDT [64]DDT [64] [64] DDT DDT 405.405 DDT [64] 405.405 [64] DDT [64] 405.405 DDT [64] 405.405 405.405 405.405 405.405 405.405 405.405 405.405 405.405

113 113 113 113 113 9.162 9.162 9.162 9.162 9.162

113 113 113 [14] [14] 113 113 [14] [14] [14] [14] 9.162 [14] [14] [14] 9.162 [14] [14] 9.162 9.162 9.162

malathion malathion malathion malathion malathion malathion malathion malathion malathion malathion

290 290 290 290 290 290 [65] [65] 2,4-D 290 290 290 290 2,4-D 2,4-D [65] 290 2,4-D [65] 2,4-D [65]2,4-D 2,4-D [65] [65] 2,4-D 23.513 2,4-D [65] [65] 23.513 2,4-D 23.513 [65] 2,4-D [65] 23.513 23.513 23.513 23.513 23.513 23.513 23.51323.513

methyl methyl methyl methyl methyl methyl methyl methyl methyl methyl parathion parathion parathion parathion parathion parathion parathion parathion parathion parathion

18 18 18 18 18 18 dicamba [67] 18 dicamba [67] [67] 18 18 18 dicamba [67] 18 dicamba [67] dicamba [67] dicamba 1.459 dicamba [67] dicamba [67] [67] [67] dicamba dicamba 1.459 1.459 dicamba [67] 1.459 1.459 1.459 1.459 1.459 1.459 1.459 1.459

75 75 75 75 75 75 75 [68] 75 [68] [68] 75 75 [68] 75 61.378 [68] [68] [68] [68] [68] [68] 61.378 61.378 61.378 61.378 61.378[68] 61.378 61.378 61.378 61.37861.378

160 160 160 160 160 [69] [69] 160 160 DEET 160 DEETDEET 160 DEET [69] 160 12.972 160 DEET [69] DEET [69]DEET DEET [69] [69] DEET [69] [69] [69] DEET 12.972 12.972 DEET [69] 12.972 12.972 12.972 12.972 12.972 12.97212.972 12.972

1.95 1.95 1.95 1.95 1.95 0.158 0.158 0.158 0.158 0.158

2 2 2222 [71] sulfoxaflor [71] 22222 sulfoxaflor [71] sulfoxaflor [71] sulfoxaflor sulfoxaflor [71] sulfoxaflor [71] sulfoxaflor 0.162 sulfoxaflor 0.162 sulfoxaflor [71] 0.162 sulfoxaflor [71] [71] sulfoxaflor [71] 0.162 0.162[71] 0.162 0.162 0.162 0.1620.162 0.162

750 750 750 750 750 750 [72] [72] 750 750 [72] 750 750 60.810 750 [72] [72] [72] [72] 60.81060.810 [72] 60.810 [72] [72] [72] 60.810 60.810 60.810 60.810 60.810 60.81060.810

naled (dibrom) naled (dibrom) naled (dibrom) naled naled (dibrom) naled (dibrom) (dibrom) naled (dibrom) naled (dibrom) naled (dibrom) nalednaled (dibrom) (dibrom) parathion parathion parathion parathion parathion parathion parathion parathion parathion parathion parathion

Cl Cl Cl Cl Cl Cl Cl Cl

O O O O O O O O Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl Cl

O O O O O O O O OHO O O OHO O OH OH OHO OH OH OHCl Cl Cl Cl Cl Cl Cl

O O O O O O O OH OH OH OH OH OH OH

639 639 639 639 639 639 [66] 639 639 [66] 639 639 [66] 639 51.810 [66] [66] [66] [66] [66] 51.810 51.810 [66] [66] [66] 51.810 51.810 51.810 51.810 51.810 51.810 51.81051.810

1.95 1.95 1.95 1.95 [70] [70] [70] 1.95 [70] 1.95 [70] [70] 0.158 [70] [70] [70] 0.158 [70] [70] 0.158 0.158 0.158 0.158

Among the pesticides (training and test set the of Among all of ofall theof listed pesticides (training and test test set compounds), the pharmacology pharmacology of Among all oflisted the listed listed pesticides (training andset testcompounds), set compounds), compounds), the pharmacology pharmacology of Among all the pesticides (training and the of Among all of the listed pesticides (training and test set compounds), the pharmacology of Among all the listed pesticides (training and set compounds), the of Among the listed pesticides (training and test set compounds), the pharmacology of Among all of ofall theof listed pesticides (training and test test set compounds), the pharmacology pharmacology of organophosphorus (OP), compounds as hAChE inhibitors, is by far the most investigated; it involves organophosphorus (OP), compounds as hAChE inhibitors, is by far the most investigated; it involves Among all of the listed pesticides (training and test set compounds), the pharmacology organophosphorus organophosphorus (OP), (OP), compounds compounds as hAChE as hAChE inhibitors, inhibitors, is by far is by the far most the investigated; most investigated; it involves it involves organophosphorus (OP), compounds as hAChE hAChE inhibitors, is by by far the investigated; most investigated; investigated; it involves involves of organophosphorus compounds as inhibitors, is far the most it organophosphorus (OP), compounds as inhibitors, is the most it organophosphorus (OP), compounds as hAChE inhibitors, is by far the most investigated; it involves organophosphorus (OP),(OP), compounds as hAChE hAChE inhibitors, is by by far far the most investigated; it involves involves the of catalytic Ser203 and the of inactive phosphyl-AChE adducts the phosphorylation phosphorylation of catalytic Ser203 andas the formation of stable stable inactive phosphyl-AChE adducts the phosphorylation phosphorylation ofcompounds catalytic Ser203 and the formation formation of stable stable inactive phosphyl-AChE adducts it the of catalytic Ser203 and the formation of inactive phosphyl-AChE adducts organophosphorus (OP), hAChE inhibitors, is by far the most investigated; the phosphorylation of catalytic Ser203 and the formation of stable inactive phosphyl-AChE adducts the of Ser203 and formation of inactive phosphyl-AChE adducts the phosphorylation of catalytic Ser203 and the formation of stable inactive phosphyl-AChE adducts the phosphorylation phosphorylation of catalytic catalytic Ser203 and the the formation of stable stable inactive phosphyl-AChE adducts [36]. Moreover, methylcarbamate (MC) insecticides like carbofuran and similar compounds, occupy [36]. Moreover, methylcarbamate (MC) insecticides like carbofuran and similar compounds, occupy [36]. Moreover, methylcarbamate (MC) insecticides like carbofuran and similar compounds, occupy [36]. Moreover, methylcarbamate (MC) insecticides like carbofuran and similar compounds, occupy [36]. Moreover, methylcarbamate (MC) insecticides like carbofuran and similar compounds, occupy involves theMoreover, phosphorylation of (MC) catalytic Ser203 and the of stable inactive phosphyl-AChE [36]. methylcarbamate insecticides like and compounds, occupy [36]. methylcarbamate (MC) insecticides likeformation carbofuran and similar compounds, occupy [36]. Moreover, Moreover, methylcarbamate (MC) insecticides like carbofuran carbofuran and similar similar compounds, occupy the catalytic triad area and act as reversible competitive inhibitors of ACh [73], elucidated by the catalytic catalytic triad area and act as reversible competitive inhibitors of ACh [73], elucidated by the the catalytic triad triad area and area act and as act reversible as reversible competitive competitive inhibitors inhibitors of ACh of ACh [73], elucidated [73], elucidated by by the[36]. catalytic triad area andas act as reversible reversible competitive inhibitors of ACh ACh [73], elucidated by the catalytic area and as competitive inhibitors of [73], elucidated adducts Moreover, methylcarbamate (MC) insecticides like carbofuran andelucidated similar compounds, the triad area act reversible competitive inhibitors of [73], by the catalytic triad area and act as reversible competitive inhibitors of ACh [73], elucidated by the catalytic catalytic triad triad area and and act asact reversible competitive inhibitors of ACh ACh [73], elucidated by by molecular docking (SB) studies on Nephotettix cincticeps AChE. Docking studies were also performed molecular docking (SB) studies on Nephotettix cincticeps AChE. Docking studies were also performed molecular molecular docking docking (SB) studies (SB) studies on Nephotettix on Nephotettix cincticeps cincticeps AChE. AChE. Docking Docking studies studies were also were performed also performed molecular docking (SB) studies on Nephotettix cincticeps AChE. Docking studies were also performed molecular docking (SB) studies on Nephotettix cincticeps AChE. Docking studies were also performed molecular docking (SB) studies on Nephotettix cincticeps AChE. Docking studies were also performed molecular docking (SB) studies on Nephotettix cincticeps AChE. Docking studies were[73], also elucidated performed by occupy the docking catalytic triad area and act as reversible competitive inhibitors of ACh molecular (SB) studies on Nephotettix cincticeps AChE. Docking studies were also performed for imidacloprid analogues [74]. for imidacloprid imidacloprid analogues [74]. [74]. for for imidacloprid analogues analogues [74]. [74]. for imidacloprid imidacloprid analogues [74]. for analogues for analogues [74]. for imidacloprid analogues for imidacloprid imidacloprid analogues [74]. [74].

Molecules 2018, 23, 2192

6 of 37

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 Molecules 2018, 23,23, 3737 Molecules 2018, 23, ofof 37 Molecules 2018, x herein study, their mechanisms of action were either confirmed or nominated 6 37 of Molecules 23, xxx of 666of pesticides, the2018, targets Molecules 2018, 23, x 6 of 37 Molecules 2018, 23, x and reversible-like (for organophosphorus pesticides only, docking performed 6 of 37 by means of reversible Therefore, toto enrich the information related toto the pharmacology ofof commonly used commercial Therefore, to enrich the information related to the pharmacology of commonly used commercial Therefore, enrich the information related the pharmacology commonly used commercial Therefore, to enrich the information related to the pharmacology of commonly used commercial to obtain non-covalent conformation that is going tomechanisms bethe transformed to the covalent one) cross-docking, Therefore, to enrichof the information related to pharmacology ofcommonly commonly used commercial pesticides, the targets of herein study, their mechanisms ofofaction action were either confirmed oror pesticides, the targets of herein study, their of action were either confirmed or pesticides, the targets ofherein hereinstudy, study, their mechanisms action wereeither either confirmed Therefore, to enrich the information related to the pharmacology of used commercial pesticides, the targets their mechanisms of were confirmed or as wellnominated as by single step molecular dynamics studies on Mus musculus and Homo sapiens AChE level. pesticides, the targets of herein study, their mechanisms of action were either confirmed or nominated by means ofofreversible reversible and reversible-like (for organophosphorus pesticides only, nominated by means of reversible and reversible-like (for organophosphorus pesticides only, nominated bymeans meansof reversible and reversible-like (for organophosphorus pesticides only, pesticides, the targets of herein study, their mechanisms oforganophosphorus action were either confirmed or by and reversible-like (for pesticides only, nominated byprediction means of reversible andof reversible-like (for organophosphorus pesticidesto only, Alongside withperformed the of free energy binding, generated conformations, were the docking performed toto obtain non-covalent conformation that going toto bebe transformed to the docking performed to obtain non-covalent conformation that isisbioactive going to be transformed to the docking performed obtain non-covalent conformation that is going transformed to the docking to obtain non-covalent conformation that is going to be transformed the 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, asnon-covalent well asas by single step molecular dynamics studies on Mus musculus covalent one) cross-docking, as well as by single step molecular dynamics studies on Mus musculus starting points for the definition of the mechanism of action ofthat 2-chloro-1,3,5-triazine-based pesticides covalent one) cross-docking, as well by single step molecular dynamics studies on Mus musculus docking performed to obtain conformation is going tostudies be transformed to the covalent one) cross-docking, as well as by single step molecular dynamics on Mus musculus covalent one) cross-docking, asAlongside well asby bywith single step molecular dynamics studies onMus Mus musculus and Homo sapiens AChE level. Alongside with the prediction ofdynamics free energy ofof binding, generated and Homo sapiens AChE level. Alongside with the prediction of free energy of binding, generated (atrazine, simazine, and propazine, respectively), as the most toxic ones. DFT quantum mechanics and Homo sapiens AChE level. Alongside with the prediction of free energy binding, generated covalent one) cross-docking, as well as single step molecular studies on musculus and Homo sapiens AChE level. the prediction of free energy of binding, generated 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 ofof bioactive conformations, were the starting points for the definition of the mechanism of action of 2-2-2bioactive conformations, were the starting points for the definition of the mechanism of action and Homo sapiens AChE level. the Alongside with the of free energy of binding, generated bioactive conformations, were starting points forprediction the definition of the mechanism offor action of 2mechanistic studies were supposed to shed further light on pesticide-AChE interactions the purpose bioactive conformations, were the starting points for the definition of the mechanism of action of 2chloro-1,3,5-triazine-based pesticides (atrazine, simazine, and propazine, respectively), as the most chloro-1,3,5-triazine-based pesticides (atrazine, simazine, and propazine, respectively), as the most chloro-1,3,5-triazine-based pesticides (atrazine, simazine, and propazine, respectively), as the most bioactive conformations, were the starting points for the definition of the mechanism of as action of 2chloro-1,3,5-triazine-based pesticides (atrazine, simazine, and propazine, respectively), the most of designing pesticides with lower acute toxicity profiles. chloro-1,3,5-triazine-based pesticides (atrazine, simazine, and propazine, respectively), asthe the most toxic ones. DFT quantum mechanics mechanistic studies were supposed toto shed further light on toxic ones. DFT quantum mechanics mechanistic studies were supposed to shed further light on toxic ones. DFT quantum mechanics mechanistic studies were supposed shed further light on chloro-1,3,5-triazine-based pesticides (atrazine, simazine, and propazine, respectively), as most toxic ones. DFT quantum mechanics mechanistic studies were supposed to shed further light on 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. pesticide-AChE interactions for the purpose of designing pesticides with lower acute toxicity profiles. pesticide-AChE interactions purpose designing pesticides with lower acute toxicity profiles. toxic ones. DFTinteractions quantum mechanics mechanistic studies were supposed to acute shed further light on pesticide-AChE forfor thethe purpose of of designing pesticides with lower toxicity profiles. 2. Results and Discussion pesticide-AChEinteractions interactionsfor forthe thepurpose purposeof ofdesigning designingpesticides pesticideswith withlower loweracute acutetoxicity toxicityprofiles. profiles. pesticide-AChE 2. Results and Discussion 2. Results and Discussion 2. Results and Discussion 2. Results and Discussion 2.1. The Conformational Analysis of Pesticides in Aqueous Solution Resultsand andDiscussion Discussion 2.2.Results 2.1. The Conformational Analysis Pesticides Aqueous Solution 2.1. The Conformational Analysis ofof Pesticides ininherein Aqueous Solution Conformational analysis (CA)of studies arein reported for simazine [9], monocrotophos [10], 2.1. The Conformational Analysis of Pesticides in Aqueous Solution 2.1. The Conformational Analysis Pesticides Aqueous Solution 2.1. The Conformational Analysis of Pesticides in Aqueous Solution dimethoate [11], and acetamiprid [12] (Table 3), as pesticides with known [9], crystal structures,[10], as well 2.1. The Conformational Analysis of Pesticides in Aqueous Solution Conformational analysis (CA) studies are herein reported for simazine [9], monocrotophos [10], Conformational analysis (CA) studies are herein reported for simazine monocrotophos Conformational analysis (CA) studies are herein reported for simazine [9], monocrotophos [10], Conformational analysis (CA) studies are herein reported for simazine [9], monocrotophos [10], as for atrazine and carbofuran, as pesticides that exert the highest acute toxicity (Table 4). To avoid Conformational analysis (CA) studies are herein reported for simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] (Table 3), asas pesticides with known crystal structures, asas well dimethoate [11], and acetamiprid [12] (Table 3), as pesticides with known crystal structures, as well dimethoate [11], and acetamiprid [12] (Table 3), pesticides with known crystal structures, well Conformational analysis (CA)[12] studies are3), herein reportedwith for simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid (Table as pesticides known crystal structures, as well dimethoate [11], and acetamiprid [12] (Table 3), as pesticides with known crystal structures, as well redundancy, applied studies to the other pesticides are reported as Supplementary Material. as for atrazine and carbofuran, as pesticides that exert the highest acute toxicity (Table 4). To avoid as for atrazine and carbofuran, as pesticides that exert the highest acute toxicity (Table 4). To avoid as for atrazine and carbofuran, as pesticides that exert the highest acute toxicity (Table 4). To avoid dimethoate [11], acetamiprid (Table that 3), asexert pesticides with acute knowntoxicity crystal(Table structures, wellThe as for atrazine and carbofuran, as [12] pesticides the highest 4). Toasavoid 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 redundancy, applied studies to the other pesticides are reported as Supplementary Material. The crystalredundancy, structures of simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] (Table applied studies to the other pesticides are reported Supplementary Material. The 3, asredundancy, for atrazineapplied and carbofuran, pesticides that exertare the highest acute toxicity (TableMaterial. 4). To avoid studies toas the other pesticides reported as as Supplementary The 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 crystal structures of simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] (Table Tablescrystal S1–S3) were used as starting conformational geometries for the CA. For pesticides without crystal structures simazine monocrotophos [10], dimethoate [11], and acetamiprid [12] (Table redundancy, applied studies to [9], the other pesticides are reported as Supplementary Material. The structures of of simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] (Table crystal structures of simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] (Table 3,3, Tables S1–S3) were used asas starting conformational geometries for the CA. For pesticides without 3, Tables S1–S3) were used as starting conformational geometries for the CA. For pesticides without Tables S1–S3) were used starting conformational geometries for the CA. For pesticides without 3, Tables S1–S3) were used as starting conformational geometries for the CA. For pesticides without crystal structures of simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] (Table experimental data (Table 4, Tables S1–S3), starting conformations were built from scratch and energy TablesS1–S3) S1–S3) were used as starting conformational geometrieswere for the CA. Forscratch pesticides without experimental data (Table Tables S1–S3), starting conformations were built from scratch and energy experimental data (Table 4,4,as Tables S1–S3), starting conformations were built from scratch and energy experimental data (Table 4, Tables S1–S3), starting conformations were built from scratch and energy 3,3, Tables were used starting conformational geometries for the CA. For pesticides without experimental (Table 4, Tables S1–S3), starting conformations built from and energy minimized. For thedata sake of clarity, pesticides were divided into four groups based on their common experimental data (Table 4,of Tables S1–S3), starting conformations were builtfrom from scratch and energy minimized. For the sake of clarity, pesticides were divided into four groups based on their common minimized. For the sake of clarity, pesticides were divided into four groups based on their common minimized. For the sake clarity, pesticides were divided into four groups based on their common experimental data (Table 4, Tables S1–S3), starting conformations were built 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 (Tables 3 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 structural patterns: 2-chloro-1,3,5-triazine derivatives: simazine, atrazine, and propazine (Table and structural patterns: 2-chloro-1,3,5-triazine derivatives: simazine, atrazine, and propazine (Table 3and and minimized. For the 2-chloro-1,3,5-triazine sake of clarity, pesticides were divided into atrazine, four groups based on their common patterns: derivatives: simazine, and propazine (Table 333and and 4, structural Table S1, respectively); amide derivatives: monocrotophos, dimethoate, and carbofuran (Tables structural patterns: 2-chloro-1,3,5-triazine derivatives: simazine, atrazine, and propazine (Table 3 and 3 Table 4, Table S1, respectively); amide derivatives: monocrotophos, dimethoate, and carbofuran Table 4, Table S1, respectively); amide derivatives: monocrotophos, dimethoate, and carbofuran Table Table respectively); amide derivatives: monocrotophos, dimethoate, and carbofuran structural patterns: 2-chloro-1,3,5-triazine derivatives: simazine, atrazine, and propazine (Table 3 and Table 4, 4, Table S1,S1, respectively); amide derivatives: monocrotophos, dimethoate, and carbofuran and 4, Table S2, respectively), carbaryl and tebufenozide (Table S2); 6-chloropyridine-3-yl)methanamine Table Table S1,4, amide derivatives: monocrotophos, dimethoate, and carbofuran carbofuran (Table 3and and Table 4,respectively); Table S2, respectively), carbaryl and tebufenozide (Table S2); 6-chloropyridine(Table and Table 4,respectively); Table S2, respectively), carbaryl and tebufenozide (Table S2); 6-chloropyridine(Table 3Table and Table 4, Table S2, respectively), carbaryl and tebufenozide (Table S2); 6-chloropyridineTable 4, S1, amide derivatives: monocrotophos, dimethoate, and (Table 334, Table Table S2, respectively), carbaryl and tebufenozide (Table S2); 6-chloropyridinederivates: acetamiprid (Table 3, Table S3) and imidacloprid (Table S3); 1-(4-chlorophenyl)-3-methylurea (Table 3 and Table 4, Table S2, respectively), carbaryl and tebufenozide (Table S2); 6-chloropyridine3-yl)methanamine derivates: acetamiprid (Table 3, Table S3) and imidacloprid (Table S3); 1-(43-yl)methanamine derivates: acetamiprid (Table 3, Table S3) and imidacloprid (Table S3); 1-(43-yl)methanamine derivates: acetamiprid (Table Table and imidacloprid (Table S3); 1-(4(Table 3 and Table 4, Table S2, acetamiprid respectively), carbaryl and tebufenozide (Table S2); 6-chloropyridine3-yl)methanamine derivates: (Table 3, 3, Table S3)S3) and imidacloprid (Table S3); 1-(43-yl)methanamine derivates: acetamiprid (Table 3, Table S3) and imidacloprid (Table S3); 1-(4condensates: diuron, monuron, and linuron (Table S3). chlorophenyl)-3-methylurea condensates: diuron, monuron, and linuron (Table S3). chlorophenyl)-3-methylurea condensates: diuron, monuron, and linuron (Table S3). chlorophenyl)-3-methylurea condensates: diuron, and linuron (Table S3). chlorophenyl)-3-methylurea diuron, monuron, and linuron (Table S3). 3-yl)methanamine derivates:condensates: acetamiprid (Table 3,monuron, Table S3) and imidacloprid (Table S3); 1-(4chlorophenyl)-3-methylureacondensates: condensates:diuron, diuron,monuron, monuron,and andlinuron linuron(Table (TableS3). S3). chlorophenyl)-3-methylurea Table 3. Training set pesticides with the known crystal structures: chemical structures, Table 3. Training set pesticides with the known crystal structures: chemical structures, Training set pesticides with the known crystal structures: chemical structures, conformational Table 3. Table3. 3.Training Trainingsetsetpesticides pesticideswith withthetheknown knowncrystal crystalstructures: structures:chemical chemical structures, Table structures, Table 3. Training set pesticides with the known crystal structures: chemical structures, conformational analysis, superposition of experimental conformation and generated global conformational analysis, superposition of the experimental conformation and generated global analysis, superposition of the experimental conformation and generated global minima using the best conformational analysis, superposition of the experimental conformation and generated global conformational analysis, of the experimental conformation generated global Table 3. Training set superposition pesticides with known crystal structures: and chemical structures, conformational analysis, superposition of the experimental conformation and generated global minima using the best performing force fields. minima using the best performing force fields. minima using the best performing force fields. minima using the best performing force fields. performing force fields. conformational analysis, superposition of the experimental conformation and generated global minima using the best performing force fields. a c minima using the best performing force fields. bb glob_minb NGMS Pesticide FF EEglob_min glob_min NGMS Pesticide FF aa a cc c Eglob_min b NGMS NGMS Pesticide FF E Pesticide FF a b a c b Pesticide FFFFa NGMSc c E(kJ/mol) E glob_minb NGMS Pesticide (kJ/mol) (kJ/mol) glob_min (kJ/mol) NGMS E glob_min Pesticide FF h h (kJ/mol) simazine NA MM3 NA simazine NA MM3 NA h h simazine (kJ/mol) NA MM3 NA simazine NA MM3 NA (kJ/mol) h simazine NA MM3 NA AMBER94 NA AMBER94 NA hh AMBER94 NA NA AMBER94 NA NA simazine MM3 NA simazine MM3 NA NA MMFF −1022.18 136 AMBER94 NA NA MMFF −1022.18 136 MMFF −1022.18 136 MMFF −1022.18 136 AMBER94 NA NA AMBER94 NA NA MMFF −1022.18 136 MMFF −1022.18 136 MMFF − 1022.18 136 MMFFs NA NA MMFFs NA NA MMFFs NA NA MMFFs NA NA MMFFs NA NA MMFFs NA NA MMFFs NA NA OPLSAA 158 OPLSAA −793.61 158 OPLSAA −−793.61 793.61 158 OPLSAA −793.61 158 OPLSAA −793.61 158 OPLSAA −793.61 158 monocrotophos MM3 NA NA monocrotophos MM3 NA NA monocrotophos MM3 NA NA monocrotophos MM3 NA NA OPLSAA −793.61 158 monocrotophos MM3 NA NA monocrotophos MM3 NA NA AMBER94 AMBER94 NA NA AMBER94 NA NA AMBER94 NA NA AMBER94 NA NA monocrotophos MM3 MMFF 241 AMBER94 NA NA MMFF −307.36 241 MMFF −307.36 241 MMFF −307.36 241 MMFF −−307.36 307.36 241 AMBER94 NA NA MMFF −307.36 241 MMFFs −310.26 218 MMFFs −310.26 218 MMFFs − 310.26 218 MMFFs −310.26 MMFFs −310.26 218 MMFF −307.36 241218 OPLSAA NA NA MMFFs −310.26 218 OPLSAA NA NA OPLSAA NA NA OPLSAA NA NA OPLSAA NA NA MMFFs −310.26 218 OPLSAA NA NA dimethoate MM3 dimethoate MM3 NA NA MM3 NA NA dimethoate MM3 NA NA dimethoate MM3 NA NA OPLSAA dimethoate AMBER94 MM3 NA NA AMBER94 NA NA AMBER94 NA NA AMBER94 NA NA AMBER94 NA NA dimethoate MM3 MMFF −−383.68 383.68 340 AMBER94 NA NA MMFF 340 MMFF −383.68 340 MMFF −383.68 340 MMFF −383.68 340 AMBER94 NA NA MMFFs −−383.68 386.71 215 MMFFs −386.71 215 MMFF 340 MMFFs −386.71 215 MMFFs −386.71 MMFFs −386.71 215 MMFF −383.68 340215 OPLSAA −−386.71 228.66 998 OPLSAA −228.66 998 MMFFs 215 OPLSAA −228.66 998 OPLSAA −228.66 OPLSAA −228.66 998 MMFFs −386.71 215998 OPLSAA −228.66 998 OPLSAA −228.66 998 acetamiprid MM3 NA NA acetamiprid MM3 NA NA acetamiprid MM3 NA NA acetamiprid MM3 NA NA acetamiprid MM3 NA NA AMBER94 AMBER94 NA NA AMBER94 NA NA AMBER94 NA NA acetamiprid MM3 MMFF −31.45 1916 AMBER94 NA NA MMFF −31.45 1916 MMFF −31.45 1916 MMFF −31.45 1916 AMBER94 NA NA MMFF −31.45 1916 MMFF −31.45 1916

dd d NF NF d NF NF NFd d NF NF d NA NA NA NA NA NA NA NA NA 2782 NA 2782 2782 2782 NA NA 2782 2782 2782 NA NA NA NA NA NA NA 2306 2306 2306 2306 2306 2306 NA NA NA NA 2306 NA NA NA NA NA NA 62 NA 62 62 62 62 NA 62 61 61 61 61 62 61 NA 61 NA NA NA NA 61 NA NA NA NA NA NA NA NA NA NA 44 NA 44 44 44 44 NA 40 40 44 40 40 44 40 15 15 40 15 15 40 15 15 15 NA NA NA NA NA NA NA NA 1414 NA 14 14 NA 14 14

Pesticide Pesticide Pesticide Pesticide e Pesticide Alignment Alignment ee e Alignment Alignment Pesticide e Alignment Alignmente e Alignment

f CAA CAA ff f CAA CAA f f gg g CAA CAA RMSD (Å) RMSD (Å) g f (Å) RMSD RMSD CAA(Å) RMSD (Å) EC/MMFF g gg EC/MMFF EC/MMFF EC/MMFF RMSD (Å)(Å) RMSD EC/MMFF 0.667 0.667 0.667 0.667 EC/MMFF EC/MMFF 0.667 0.667 0.667

EC/MMFF EC/MMFF EC/MMFF EC/MMFF EC/MMFF EC/MMFF 0.620 0.620 0.620 0.620 0.620 EC/MMFF EC/MMFFs 0.620 EC/MMFFs EC/MMFFs EC/MMFFs EC/MMFFs 0.620 EC/MMFFs 3.521 3.521 3.521 3.521 3.521 EC/MMFFs 3.521 3.521 EC/MMFF EC/MMFF EC/MMFF EC/MMFF EC/MMFF 0.148 EC/MMFF 0.148 0.148 0.148 0.148 EC/MMFF EC/MMFFs 0.148 EC/MMFFs EC/MMFFs EC/MMFFs EC/MMFFs 0.148 3.227 3.227 EC/MMFFs 3.227 3.227 3.227 EC/MMFFs EC/OPLSAA EC/OPLSAA 3.227 EC/OPLSAA EC/OPLSAA EC/OPLSAA 3.227 2.746 EC/OPLSAA 2.746 2.746 2.746 2.746 EC/OPLSAA 2.746 EC/MMFF EC/MMFF EC/MMFF EC/MMFF 2.746 EC/MMFF 0.783 0.783 0.783 0.783 EC/MMFF EC/MMFFs 0.783 EC/MMFFs EC/MMFFs EC/MMFFs 0.783 EC/MMFFs EC/MMFFs

Molecules 2018, 23, 2192

7 of 37

Table 3. Cont. Pesticide

Molecules 2018, 23, x Molecules 2018, Molecules 2018, 23,23, x x Molecules 2018, Molecules 2018, Molecules 2018, 23,23, x23,x x

FF a

Eglob_min (kJ/mol)

b

NGMS c

NF d

Pesticide Alignment e

CAA f

7 of of 37 37 7 of7 37

g ofof3737 RMSD (Å) 7 of7 737

MM3 NA NA NA EC/MMFF AMBER94 NA NA NA 0.783 MMFFs −31.29 −31.29 1816 6 0.823 MMFFs −31.29 1816 0.823 MMFFs 0.823 MMFF −31.45 1816 1916 6 6 14 EC/MMFFs MMFFs −31.29 1816 6 0.823 MMFFs −31.29 1816 6 0.823 MMFFs −31.29 1816 6 0.823 MMFFs −31.29 1816 6 0.823 OPLSAA NA NA NA OPLSAA NA NA NA OPLSAA NA c NA NA NA NA NA a Force field; b Energy OPLSAA the global minimum; Number of that aofsingle global minimum structure was found OPLSAA NA NA NA a Force field; of b Energy ctimes OPLSAA NAminimum; NA NA of the global Number times that a single global minimum OPLSAA NA NA NA a Force b Energy c Number a Force b Energy c Number field; of the global minimum; of times that a single global minimum field; of the global minimum; of times that a single global minimum d in 10,000 processed structures; of families i.e.,c different conformations found in 10,000 processed structures; a aForce b bEnergy ofNumber c cNumber d Number field; the global minimum; of times that a single global minimum Force field; Energy of the global minimum; Number of times that a single global minimum a Force b of families i.e., different conformations structure was found in 10,000 processed structures; d Number field; Energy of 10,000 the processed global minimum; Number ofBlue—MMFFs thati.e., ai.e., single global minimum d Number e Red—experimental of families different conformations structure was found processed structures; oftimes families different conformations structure was found in in 10,000 structures; conformation, Violet—MMFF conformation, conformation, Green—OPLSAA d d e dRed—experimental Number ofoffamilies i.e., different conformations structure found inin 10,000 processed structures; Number families i.e., different conformations structure was found 10,000 processed structures; found in 10,000 processed structures; conformation, Violet—MMFF f was g RMSD Number of families i.e., different conformations structure found in processed 10,000 processed structures; e eRed—experimental found in 10,000 structures; Red—experimental conformation, Violet—MMFF conformation; Conformation dependent alignment accuracy; measured between the heavy atoms of found inwas 10,000 processed structures; conformation, Violet—MMFF e e Red—experimental fViolet—MMFF found inin 10,000 processed structures; conformation, Violet—MMFF found 10,000 processed structures; Red—experimental conformation, e Red—experimental h Conformation conformation, Blue—MMFFs conformation, Green—OPLSAA conformation; f found in 10,000 processed structures; conformation, Violet—MMFF pesticides of experimental and best performing force Green—OPLSAA field conformations;conformation; Not available.f Conformation Conformation conformation, Blue—MMFFs conformation, Green—OPLSAA conformation; conformation, Blue—MMFFs conformation, f f Conformation g RMSD measured conformation, Blue—MMFFs conformation, Green—OPLSAA conformation; Conformation conformation, Blue—MMFFs conformation, Green—OPLSAA conformation; f between the heavy atoms ofpesticides pesticides dependent alignment accuracy; g RMSD conformation, Blue—MMFFs conformation, Green—OPLSAA conformation; g RMSD measured between heavy atoms of dependent alignment accuracy; measured between thethe heavy atoms of Conformation pesticides of of of dependent alignment accuracy; g g h RMSD measured between the heavy atoms of dependent alignment accuracy; RMSD measured between the heavy atoms of pesticides pesticides dependent alignment accuracy; g RMSD experimental and best performing force field conformations; Not available. h measured between the heavy atoms of pesticides of ofof dependent alignment accuracy; h and best performing force field conformations;Not Not available. best performing force field conformations; available. hpesticide CAexperimental inexperimental explicit and solvent may be used to predict the single stability prior to AChE binding. h Notavailable. experimental and best performing force field conformations; Not experimental and best performing force field conformations; available. h Not experimental and best performing force field conformations; available. acetamiprid

Herein, CAs were tentatively carried out by using five different force fieldsprior (FFs): MM3, AMBER94, CA in explicit solvent may be used to predict the single pesticide stability to AChE binding. CA in explicit solvent may used predict single pesticide stability AChE binding. CA in explicit solvent may bebe used to to predict thethe single pesticide stability priorprior to to AChE binding. CA in explicit solvent may be used to predict the single pesticide stability prior to AChE binding. 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 CA in explicit solvent may be used to predict the single pesticide stability prior to AChE binding. Herein, CAs were tentatively carried using five different force fields (FFs): MM3, AMBER94, Herein, CAs were tentatively carried outout byby using five different force fields (FFs): MM3, AMBER94, Herein, CAs were tentatively carried out by five different force fields (FFs): MM3, AMBER94, Herein, CAs were tentatively carried out byusing using five different force fields (FFs): MM3, AMBER94, MMFF, MMFFs, and OPLSAA [35]. Either for crystallized pesticides or pesticides with unknown Herein, CAs were tentatively carried out by using five different force fields (FFs): MM3, AMBER94, experimental structures, the procedure was divided into the following steps: (1) conformational MMFF, MMFFs, and OPLSAA [35]. Either crystallized pesticides pesticides with unknown MMFF, MMFFs, and OPLSAA [35]. Either forfor crystallized pesticides oror pesticides with unknown MMFF, MMFFs, and OPLSAA [35]. Either for crystallized pesticides pesticides with unknown MMFF, MMFFs, and OPLSAA [35]. Either for crystallized pesticides orand pesticides with unknown experimental structures, theprocedure procedure was divided into thefollowing following steps: (1)conformational conformational MMFF, and OPLSAA [35]. Either for crystallized pesticides or or pesticides with unknown analysis; (2)MMFFs, determination ofthe energy of the global minimum (E (3) the LB superposition experimental structures, was divided into the steps: experimental structures, the procedure was divided into the following (1)(1) conformational glob_min ),steps: experimental structures, the procedure was divided into the following steps: (1) conformational experimental structures, the procedure was divided into the following steps: (1) conformational analysis; (2) determination of energy of the global minimum (E glob_min),steps: and (3) the LB superposition experimental structures, the procedure was divided into the following (1) conformational analysis; (2) determination energy of the global minimum (E glob_min ), and the LB superposition analysis; (2)to determination of of energy of the global minimum (Eglob_min ), and (3)(3) the LB superposition of structures define the alignment accuracy. All steps were performed in order to obtain the best analysis; (2) of energy of the global minimum (E(Eglob_min ),),and (3) the LB superposition analysis; (2)determination determination of energy of the global minimum glob_min and (3)order the LB superposition of structures to define the alignment accuracy. All steps were performed in to obtain the best analysis; (2) determination of energy of the global minimum (E glob_min ), and (3) the LB superposition structures to define alignment accuracy. All steps were performed order obtain best of of structures to define thethe alignment accuracy. All steps were performed in in order to to obtain thethe best performing FF, either by means of capability to predict the free energy of global minimum or to of structures to the alignment accuracy. All steps were performed in order to the best of structures todefine define the alignment accuracy. All steps were performed inof order toobtain obtain the best performing FF, either by means of capability to predict the free energy global minimum or to 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 performing FF, either by means of capability to predict the free energy of global minimum or to superimpose molecules. The goal of CA was to answer the question: tominimum consider or the FF of performing FF, either by means of to predict the free energy of global performing FF, either byhigher means ofcapability capability towas predict the free energy ofHow global minimum ortoto superimpose molecules. The higher goal of CA to answer the question: How to consider the FF performing FF, either by means of capability to predict the free energy of global minimum or to superimpose molecules. The higher goal CA was answer the question: How to consider superimpose molecules. The higher goal of of CA was to to answer theor question: How tocapability? consider thethe FFFF choice while treating pesticides: by predicted ECA alignment superimpose molecules. The higher goal of was to answer question: How to glob_min superimpose molecules. The higher goal ofCA was tovalues answer the question: How toconsider consider the FF of choice while treating pesticides: by Eglob_min values orby by alignment superimpose molecules. The higher goal ofpredicted CA was to answer thethe question: How tocapability? consider thethe FFFF choice while treating pesticides: predicted Eglob_min values alignment capability? of of choice while treating pesticides: byby predicted Eglob_min values oror byby alignment capability? while treating pesticides: EEglob_min values capability? ofchoice choice while treating pesticides: bypredicted predicted glob_min values orby byalignment alignment capability? of of choice while treating pesticides: byby predicted Eglob_min values or or by alignment capability? Table 4. Training set pesticides of the highest toxicity with the unknown crystal structures: chemical Table 4. Training set pesticides of the highest toxicity with the unknown crystal structures: chemical Table 4. Training set pesticides of the highest toxicity with the unknown crystal structures: chemical Table 4. Training set pesticides of the highest toxicity with the unknown crystal structures: chemical Table 4.4.Training set of toxicity with the crystal structures: chemical Table Training setpesticides pesticides ofthe thehighest highest toxicity with theunknown unknown crystal structures: chemical structures, conformational analysis, superposition of the experimental conformation and generated structures, conformational analysis, superposition of the experimental conformation and generated 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 structures, conformational analysis, superposition of the experimental conformation and generated structures, conformational analysis, superposition of experimental conformation and generated structures, conformational analysis, superposition ofthe the experimental conformation and generated global minima using the best performing force fields. structures, conformational analysis, superposition of the experimental conformation and generated global minima using the best performing force fields. global minima using best performing force fields. global minima using thethe best performing force fields. global minima using force global minima using thebest bestperforming performing forcefields. fields. global minima using thethe best force a performing c d b fields. Eglob_min NGMS Pesticide Pesticide a d a FF c c NF d NF b b NGMS Eglob_min NGMS NF Pesticide Pesticide E glob_min Pesticide Pesticide FFaFF cc c b dd d aa bb Pesticide FF NGMS Pesticide e E NF E glob_min NGMS NF Pesticide Pesticide FF E glob_min NGMS NF Pesticide Pesticide FF glob_min a c d b (kJ/mol) Alignment E glob_min NGMS NF Pesticide Pesticide FF (kJ/mol) Alignment (kJ/mol) Alignment e ee e e (kJ/mol) Alignment (kJ/mol) Alignment e atrazine MM3 −1161.48 445 23 (kJ/mol) Alignment (kJ/mol) Alignment atrazine MM3 −1161.48 atrazine MM3 −1161.48 445445 23 23 h atrazine MM3 −1161.48 445 2323 atrazine MM3 −1161.48 445 AMBER94 NA NA NA h hNA atrazine MM3 −1161.48 445 23NA AMBER94 NA NA AMBER94 NA atrazine MM3 −NA 1161.48 445 23 hh AMBER94 NA NA NA AMBER94 NA NA NA hh MMFF −1007.32 −1007.32 147 NA 62 AMBER94 NA NA MMFF −1007.32 147 62 MMFF 147 62NA AMBER94 NA NA MMFF −1007.32 147 6262 MMFF −1007.32 147 MMFFs −992.63 509 MMFF −1007.32 147 62 MMFF − 1007.32 147 62 MMFFs −992.63 MMFFs −992.63 509509 20 2020 MMFFs −992.63 509 2020 MMFFs −992.63 −992.63 509 MMFFs 509 2020 MMFFs −992.63 509 OPLSAA −783.57 364 OPLSAA −783.57 OPLSAA −783.57 364364 19 1919 OPLSAA − 783.57 364 19 OPLSAA −783.57 OPLSAA −783.57 −783.57 364 OPLSAA 364364 19 1919

f FFDAA f f FFDAA FFDAA ff f g FFDAA FFDAA FFDAA f RMSD g g (Å) FFDAA RMSD RMSD (Å)(Å) g RMSD (Å) RMSD (Å) g g g MM3/MMFF RMSD (Å) RMSD (Å) MM3/MMFF MM3/MMFF MM3/MMFF MM3/MMFF 0.892 MM3/MMFF 0.892 0.892 MM3/MMFF 0.892 0.892 MM3/MMFFs 0.892 MM3/MMFFs MM3/MMFFs 0.892 MM3/MMFFs MM3/MMFFs 0.931 MM3/MMFFs MM3/MMFFs 0.931 0.931 0.931 0.931 0.931 0.931 S MMFF/MMFF MMFF/MMFF S S MMFF/MMFF MMFF/MMFF SS S MMFF/MMFF MMFF/MMFF MMFF/MMFF 0.343S 0.343 0.343 0.343 0.343 0.343 carbofuran S MM3 29.60 3617 2 MMFF/MMFF 0.343 carbofuran MM3 29.60 3617 MMFF/MMFF carbofuran S S MM3 29.60 3617 2 22 MMFF/MMFF carbofuran MM3 29.60 3617 MMFF/MMFF SS carbofuran MM3 29.60 3617 2NA MMFF/MMFF carbofuran S MM3 29.60 3617 2 MMFF/MMFF AMBER94 NA NA 0.075 carbofuran S MM3 29.60 3617 2 MMFF/MMFF AMBER94 NA NA NA 0.075 AMBER94 NA NA NA 0.075 AMBER94 NA NA NA 0.075 AMBER94 NA NA NA 0.075 AMBER94 NA NA NA 0.075 MMFF 20.30 3307 NA MMFF/OPLSAA AMBER94 NA NA 0.075 MMFF 20.30 3307 MMFF/OPLSAA MMFF 20.30 3307 MMFF/OPLSAA MMFF 20.30 3307 2 22 2 MMFF/OPLSAA MMFF 20.30 3307 2 MMFF/OPLSAA MMFF 20.30 3307 2 MMFF/OPLSAA MMFFs 15.92 3384 2 0.216 MMFF 20.30 3307 MMFF/OPLSAA MMFFs 15.92 3384 0.216 MMFFs 15.92 3384 0.216 MMFFs 15.92 3384 22 22 0.216 MMFFs 15.92 3384 2 275 0.216 MMFFs 15.92 3384 0.216 OPLSAA 73.27 5621 MMFFs/OPLSAA OPLSAA 73.27 5621 75 MMFFs/OPLSAA MMFFs 15.92 3384 2 0.216 OPLSAA 73.27 5621 MMFFs/OPLSAA OPLSAA 73.27 5621 75 75 MMFFs/OPLSAA OPLSAA 73.27 5621 MMFFs/OPLSAA OPLSAA 73.27 73.27 5621 MMFFs/OPLSAA 0.239 0.239 OPLSAA 5621 75 7575 MMFFs/OPLSAA 0.239 0.239 0.239 a Force b Energy of the globalc minimum; c Number of times that a single global minimum a Force field; b Energy 0.239 field; 0.239 of the global minimum; Number of times that a single global minimum structure was a b c a Force b Energy c Number Force field; Energy global minimum; Number times that a single global minimumfound field; of of thethe global minimum; of of times that a single global minimum a aForce field; b bEnergy d c cNumber d Number of the global minimum; of that ai.e., global minimum Force field; Energy of the global minimum; Number oftimes times that asingle single global minimum a Force bwas c different in 10,000 processed structures; Number of families i.e., conformations found in 10,000 processed structures; of families different conformations structure found in 10,000 processed structures; d field; Energy of the global minimum; Number of times that a single global minimum d Number families different conformations structure was found 10,000 processed structures;Number of of families i.e.,i.e., different conformations structure was found in in 10,000 processed structures; e Red—experimental d Number d Number conformation, Violet—MMFF conformation, Blue—MMFFs conformation, Green—OPLSAA e dRed—experimental ofoffamilies i.e., different conformations structure was found inin 10,000 processed structures; families i.e., different conformations structure was found 10,000 processed structures; found in 10,000 processed structures; conformation, Violet—MMFF e Number of families i.e., different conformations structure was found in 10,000 processed structures; e foundin fin10,000 10,000processed processedstructures; structures; Red—experimental Red—experimental conformation,Violet—MMFF Violet—MMFF found conformation, g RMSD measured e e Red—experimental conformation; Force field dependent alignment accuracy; between thef Force heavy atoms of pesticides found 10,000 processed structures; conformation, Violet—MMFF found 10,000 processed structures; Red—experimental conformation, Violet—MMFF e Red—experimental field dependent conformation, Blue—MMFFs conformation, Green—OPLSAA conformation; f Force found in inin 10,000 processed structures; conformation, Violet—MMFF f Force field dependent conformation, Blue—MMFFs conformation, Green—OPLSAA conformation; h field dependent conformation, Blue—MMFFs conformation, Green—OPLSAA conformation; experimental and bestBlue—MMFFs performing force field conformations; Not available. f Force f g field dependent conformation, conformation, Green—OPLSAA conformation; Force dependent conformation, Blue—MMFFs conformation, Green—OPLSAA conformation; f Force RMSD measured between the heavy atoms of pesticides experimental and best alignment accuracy; g RMSD fieldfield dependent conformation, Blue—MMFFs conformation, Green—OPLSAA conformation; g RMSD measured between heavy atoms pesticides experimental and best alignment accuracy; measured between thethe heavy atoms of of pesticides experimental and best alignment accuracy; g RMSD g RMSDmeasured hbetween between the atoms ofofpesticides experimental and best alignment accuracy; measured theheavy heavy atoms pesticides experimental and best alignment accuracy; g RMSD performing force field conformations; Not available. h measured between the heavy atoms of pesticides experimental and best alignment accuracy; h performing force field conformations;Not Not available. performing force field conformations; available. h Not h Notavailable. performing force field conformations; performing force field conformations; available. h Not CAperforming results for simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] are force field conformations; available.

CA results forsimazine simazine [9],monocrotophos monocrotophos [10], dimethoate and acetamiprid [12] reportedCA in Table 3. MM3 and[9], AMBER94 could not be applied to all[11], ofand the crystallized pesticides CA results [9], [10], dimethoate [11], acetamiprid [12] results forfor simazine monocrotophos [10], dimethoate [11], and acetamiprid [12] areareare CA results for simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] are 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 CA results for simazine [9], monocrotophos [10], dimethoate [11], and acetamiprid [12] are reported Table MM3 and AMBER94 could applied of crystallized pesticides reported in in Table 3. 3. MM3 and AMBER94 could notnot bebe applied to to allall of thethe crystallized pesticides as as as the needed parameters for the optimization were not available. Possible workaround could reported in Table 3. MM3 and AMBER94 could not be applied to all of the crystallized pesticides as 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 reported in Table MM3 and AMBER94 notnot applied allPossible of theworkaround crystallized pesticides asthe the needed parameters optimization were not available. workaround could the needed parameters forfor thethe optimization were available. Possible could bebe the be the utilization of3. Antechamber suite could [75], by be which thetogeneric AMBER parameters could be the needed parameters the optimization were not available. Possible workaround could the needed parameters for the optimization were not available. Possible workaround could bethe the utilization of Antechamber suite [75], by which the generic AMBER parameters could be calculated, theutilization needed parameters forfor thesuite optimization were not available. Possible workaround could bebe the of Antechamber [75], by which the generic AMBER parameters could be calculated, utilization of Antechamber suite [75], by which the generic AMBER parameters could be calculated, calculated, but then a different protocol would have been developed instead of the one based on the utilization Antechamber suite [75], by the generic AMBER parameters could be utilization of Antechamber suite [75], bywhich which the generic AMBER parameters could becalculated, calculated, but then aof different protocol would been developed instead of the one based on the well-known utilization of Antechamber suite [75], byhave which the generic AMBER parameters could be calculated, but then a different protocol would have been developed instead the one based on the well-known but then a accredited different protocol would have been developed instead of of the one based on the well-known well-known Schrödinger software [76]. MMFF showed to be the only FF containing but then a different protocol would have been developed instead of the one based on the well-known but then a different protocol would have been developed instead of the one based on the well-known accredited Schrödinger software [76]. developed MMFF showed tothe beone the only containing thethe but then a different protocol would have been instead of based on FF the well-known accreditedSchrödinger Schrödingersoftware software[76]. [76].MMFF MMFFshowed showedto tobebethetheonly onlyFFFFcontaining containingthethe accredited accredited Schrödinger software [76]. MMFF showed only containing accredited Schrödinger software [76]. MMFF showed be the only FF containing the accredited Schrödinger software [76]. MMFF showed to toto bebe thethe only FFFF containing thethe

Molecules 2018, 23, 2192

8 of 37

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 (Tables 3 and 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 (Tables 3 and 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, Tables 3 and 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

Molecules 2018, 23, 2192

9 of 37

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)): pLD50 = −0.0021Eglob_min − 2.82 (r2 = 0.78),

(2)

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)): pLD50 = −0.0019Eglob_min − 3.05 (r2 = 0.905), (3) 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 (Tables 2 and 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.

Molecules 2018, 23, 2192

10 of 37

Molecules 2018, 23, x

10 of 37

Figure 2. The plot of experimental vs. fitted pLD50 values calculated with QSAR model 2 for the 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 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 indicated by orange squares) (a); test set pesticides (green circles) (b) against Mus musculus. The plot of values calculated with the QSAR model 4 for the training set pesticides of calculated calculated vs. vs. fitted fitted pLD pLD50 50 values calculated with the QSAR model 4 for the training set pesticides (training (training set set values values are are depicted depicted by by blue blue circles; circles; internal internal test test set set values values are are indicated indicated by by orange orange squares squares (c); test set pesticides (green circles) (d) against Homo sapiens. (c); test set pesticides (green circles) (d) against Homo sapiens.

Moreover, models 1 and 2 cannot be universally applied for the prediction of the Homo sapiens 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 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 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 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 set, pesticides did not lead to any statistically significant model. As above, exclusion of carbofuran and linuron led to the model 3 (Equation (4)): and linuron led to the model 3 (Equation (4)): = −0.0021 + 3.01 = 0.783 , (4) _ pLD50 = −0.0021Eglob_min + 3.01 (r2 = 0.783), (4) and further exclusion of diuron and monuron led to model 4 (Equation (5), Tables S7 and S8, Figure 2c,d): and further exclusion of diuron and monuron led to model 4 (Equation (5), Tables S7 and S8, = −0.0019 + 4.15 = 0.91 , Figure 2c,d): (5) _ pLD50 = −0.0019Eglob_min + 4.15 (r2 = 0.91), (5) Like for Mus musculus and Homo sapiens, the pLD50 values were not fully satisfactory correlated to the global minima energies, leading to the thepLD conclusion glob_min are not adequate Like for Mus musculus and Homo sapiens, wereEnot fullyvalues satisfactory correlated to 50 values that indicators to anticipate the pLD 50 against model organism, the level fromindicators LB to SB, the global minima energies, leading to theany conclusion that Eglob_min valueswas arechanged not adequate directing thatthethe free energies binding could the be level considered as appropriate toxicity to anticipate pLD anyofmodel organism, was changed from LB toacute SB, directing 50 against indicators. energies of binding be calculated separately for each pesticide and model that the freeFree energies of binding couldcan be considered as appropriate acute toxicity indicators. Free organism, performing studies, either for Mus AChE orby Homo sapiens energies ofby binding can be molecular calculated docking separately for each pesticide andmusculus model organism, performing AChE. 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

Molecules 2018, 23, 2192

11 of 37

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

Molecules 2018, 23, 2192

12 of 37

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 Molecules 2018, 23, xof structurally different inhibitors. Within the RCCD stage, Vina was far superior 12 of 37 to the presence in comparison with other programs, by means of The SB alignment accuracy, and, therefore, can be considered as as aa valuable valuable tool tool for for the the prediction prediction of of bioactive bioactive conformations conformations of of acetylcholine acetylcholine esterase esterase considered inhibitors with the unknown binding mode. inhibitors with the unknown binding mode.

2.5.2. The Cross-Docking and Molecular Molecular Dynamics Dynamics Studies Studies of of ACh ACh Both the training set or the test set pesticides were submitted to the herein Vina-based proposed been validated by reproducing the bioactive conformations of PDBSB alignment alignment protocols, protocols,that thathave have been validated by reproducing the bioactive conformations of available inhibitors in complex with AChE (Supplementary Materials Tables Tables S9 andS9 S10, and Figure PDB-available inhibitors in complex with AChE (Supplementary Materials and S10, and S2) andS2) applied in situ in situ predicting the bioactive conformations of AChof(Figure 3). Figure and applied in predicting the bioactive conformations ACh (Figure 3). no experimental experimental structures of ACh ACh in in complex complex with with either either wild-type wild-type mAChE or hAChE are As no available in the Protein Protein Data Bank, its bioactive conformation within both enzymes was predicted by available means of Vina, Vina, as as the the best performing SB alignment assessment tool (Figure 3a,b, see Supplementary Materials Alignment AlignmentRules RulesAssessment Assessment section, Tables S9 S10). and S10). The stability ofACh-mAChE either AChMaterials section, Tables S9 and The stability of either mAChE or ACh-hAChE complex was confirmed using one step molecular dynamics (MD) runs, upon upon or ACh-hAChE complex was confirmed using one step molecular dynamics (MD) runs, inspection of Cα Cα RMSD RMSD values values (1.2 (1.2 ns ns simulation simulation length, length, Figures Figures S7a S7a and and S8a), S8a), RMSF RMSF values values the inspection (Figures S7b S7b and and S8b), S8b), and and trajectories trajectories energy energy terms terms (data (data not not shown shown but but available available from from the the authors authors (Figures upon request). request). The short-termed MD protocol was sufficient to to produce produce stable stable complexes. complexes. The The two two ACh poses poses (in (in mAChE mAChE or or hAChE) hAChE) were were conformationally conformationally very very close, close, displaying displaying a RMSD RMSD value value of of ACh Å. 0.043 Å.

Figure Figure 3. 3. The SB alignment alignment of of acetylcholine acetylcholine into into mAChE mAChE active active site site (a); (a); hAChE hAChE active active site site (b). (b). The The enzyme ribbons are presented in blue and gold, respectively, respectively, active active site site amino amino acids acids are are depicted depicted in in white. For the thepurpose purposeofofclarity, clarity,hydrogen hydrogenatoms atoms omitted from presentation. graphic white. For areare omitted from thethe presentation. The The graphic was was generated using UCSF Chimera software(Resource (Resourcefor for Biocomputing, Biocomputing, Visualization, Visualization, and generated using the the UCSF Chimera software and Informatics (RBVI), University of California, San Francisco, CA, USA). Informatics (RBVI), University of California, San Francisco, CA, USA).

expected, the positive quaternary ACh amine group was located in the anionic sub-site, As expected, where the orientation of of particular particular group group is is supported supported by by strong strong hydrophobic hydrophobic interactions interactions with a where the orientation m(h)Trp86 methyl methyl group, group, as as well as with m(h)Tyr124, m(h)Tyr337, m(h)Tyr337, and and m(h)Tyr341 m(h)Tyr341 aromatic aromatic moieties moieties m(h)Trp86 (Figures S7c S7c and and S8c). S8c). The The positively positively charged charged amine amine was was electrostatically electrostatically attracted attracted by the the electronic electronic (Figures cloud of the m(h)Trp86 m(h)Trp86 aromatic aromatic group. group. The The acetyl acetyl group group occupied occupied the the esterase esterase sub-site: sub-site: the carbonyl cloud interactions with m(h)Tyr337 sideside chain hydroxyl group and oxygen was was involved involvedininthe theelectrostatic electrostatic interactions with m(h)Tyr337 chain hydroxyl group m(h)His447 imidazole ring,ring, while thethe carbonyl oxygen was of and m(h)His447 imidazole while carbonyl oxygen wasoriented orientedtowards towardsthe the side side chain of moiety was observed to m(h)Tyr133. Regarding Regarding the the mouse mouseenzyme, enzyme,during duringthe theMD MDrun, run,the thecarbonyl carbonyl moiety was observed 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

Molecules 2018, 23, 2192

13 of 37

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 (Figures 4 and 5, respectively). The graphical interpretation of the results of the SB and MD studies are herein shown for atrazine (Figures 4a and 5a, Figure S9 and Figure S10, respectively) and simazine (Figures 4b and 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).

Molecules 2018, 23, 2192 Molecules 2018, 23, x

14 of 37 14 of 37

Figure SBSB alignment of atrazine (a); simazine (b); carbofuran (c); monocrotophos (d); dimethoate Figure4.4.The The alignment of atrazine (a); simazine (b); carbofuran (c); monocrotophos (d); (e); imidacloprid (f) into the mAChE site. The enzyme ribbons are presented blue, the in active dimethoate (e); imidacloprid (f) into active the mAChE active site. The enzyme ribbons areinpresented site amino acidssite are amino depicted in are white. For the tothe clarify, the hydrogen atoms are omitted. blue, the active acids depicted in purpose white. For purpose to clarify, the hydrogen atoms The are omitted. The graphic generated using the UCSF Chimera software (Resource for graphic was generated usingwas the UCSF Chimera software (Resource for Biocomputing, Visualization, Biocomputing, and Informatics (RBVI), University CA, of California, San Francisco, CA, and InformaticsVisualization, (RBVI), University of California, San Francisco, USA). USA).

Molecules 2018, 23, 2192 Molecules 2018, 23, x

15 of 37

15 of 37

Figure 5. The SB alignment of atrazine (a); simazine (b); carbofuran (c); monocrotophos (d); dimethoate Figure 5. The SB alignment of atrazine (a); simazine (b); carbofuran (c); monocrotophos (d); (e);dimethoate imidacloprid into the hAChE active site. The enzyme are presented inpresented gold; the active (e);(f) imidacloprid (f) into the hAChE active site. ribbons The enzyme ribbons are in gold;site amino acids are in white. For theinpurpose to clarify, the hydrogen are omitted from the active site depicted amino acids are depicted white. For the purpose to clarify,atoms the hydrogen atoms arethe presentation. The the graphic was generated using thewas UCSF Chimerausing software for Biocomputing, omitted from presentation. The graphic generated the (Resource UCSF Chimera software Visualization, andBiocomputing, Informatics (RBVI), University California,(RBVI), San Francisco, CA,of USA). (Resource for Visualization, andofInformatics University California, San Francisco, CA, USA).

Nevertheless, within either mice or human bioactive conformations, the 2-chloro-1,3,5-triazines Nevertheless, within either mice or human bioactive conformations, 2-chloro-1,3,5-triazines were parallel to m(h)Trp86 and m(h)Tyr337, trapped in the hydrophobic cagethe between two residues. The were parallel to m(h)Trp86 and m(h)Tyr337, trapped in the hydrophobic cage between two residues. individual contribution of these interactions was respectively −37.34 and −19.9 kcal/mol, according The individual contribution of these interactions was respectively −37.34 and −19.9 kcal/mol, to the Free-Energy Decomposition Analysis, FEDA (Table S11). Atrazine bioactive conformations according to the Free-Energy Decomposition Analysis, FEDA (Table S11). Atrazine bioactive were nearly identical in both enzymes. On the other hand, comparing simazine docked structures (Figures 4b and 5b), the one found in hAChE (Figure 5b) was in-plane rotated for about 30◦ . The

Molecules 2018, 23, 2192

16 of 37

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 0 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 0 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 (Figures 4c and 5c, Figure S15 and Figure S16, respectively) and monocrotophos (Figures 4d and 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 (Figures 4c and 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

Molecules 2018, 23, 2192

17 of 37

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 (Figures 4e and 5e, respectively). Hence, although dimethoate is in both active sites superimposed to monocrotophos, by means of backbone orientation (Figures 4d and 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.

Molecules 2018, 23, 2192

18 of 37

Molecules 2018, 23, x

21 of 37

Figure 6. Schematic view of pesticides binding interactions within the hAChE active site: atrazine, Figure 6.or Schematic of pesticides within the hAChE active atrazine, simazine propazineview (a); carbofuran or binding carbarylinteractions (b); monocrotophos or dimethoate (c) site: imidacloprid 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 or acetamiprid (d); diuron, monuron,interactions or linuron (e). hydrogen arethe presented with pink dashed half arrows, the hydrophobic withThe a blue dot-linebonds marker, ionic interactions 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 with thehAChE red dashed against (f). lines. The proposed mechanism of monocrotophos or dimethoate acute toxicity against hAChE (f).

Carbaryl reached the comparable bioactive conformation to carbofuran (Figures S3b and S5b), The results of molecular docking are presented as a Supplementary Materials (Figures S35–S42). within both species AChE active sites, with the difference that the 2,2-dimethyltetrahydrofuran scaffold The generated conformations were used further on for the external validation of the training set is in the structure of particular pesticide replaced with another benzene ring. Such replacement has a results.

Molecules 2018, 23, 2192

19 of 37

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-3Yl)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 (Figures 4f and 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 (Figures 4f and 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).

Molecules 2018, 23, 2192

20 of 37

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 0 -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 diuron0 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

Molecules 2018, 23, 2192

21 of 37

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): pLD50 = −0.7885∆Gbinding − 2.44 (r2 = 0.95),

(6)

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.



= −0.7885

binding

− 2.44

= 0.95 ,

(6)

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 humans on the basis of ∆Gbinding values for Homo sapiens AChE (Figure227c,d, Molecules 2018, 23,against 2192 of 37 Tables S14 and S15).

Figure 7. The plot of experimental vs. fitted pLD50 values calculated with QSAR model 5 for training 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 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 values calculated with QSAR model 5 for training set pesticides (c); test set pesticides (d) against Homo sapiens. sapiens.

Thus, the QSAR model 5 described the acute toxicities against Mus musculus, by means of the 2.7. The Quantum Mechanical Studies of Acetylcholinesterase Inhibition by 2-Chloro-1,3,5-triazine-based free energies of binding of reversible and reversible-like bioactive conformations for either nonPesticides covalent or covalent pesticides, respectively. The model is also proven to be a tool for the anticipation Thetoxicities performed molecular andthemolecular dynamics studies agreed with the of acute against humans.docking In the end, conduced structure-based studies unequivocally previously reported mode of action of amide-based, 6-chloropyridine-3-yl)methanamine-based, supported pLD50/∆Gbinding interrelation obtained by QSAR model 5. In that sense, the SB approach, and 1-(4-chlorophenyl)-3-methylurea-based pesticides The quantum mechanical (QM) conducted through Vina as the cheminformatics engine, [73]. was proven to be the correct one for the calculations were only performed on hAChE, with the aim to reveal the pharmacology of anticipation of acute toxicities of targeted pesticides. Consequently, the Vina-based ∆Gbinding values 2-chloro-1,3,5-triazine-based pesticides and theindicators origin for acute the toxicity, givenLD that their can be considered in future design as valuable to the anticipate pesticides 50 values pharmacology is the least understood. against any model organism. 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)

Molecules 2018, 23, 2192

23 of 37

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-N0 atom hydrogen towards hThr83 (Figure 8a). The distance between the hThr83 hydroxyl group oxygen atom as a HB acceptor and m-N0 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-N0 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-N0 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 (Figures 8a and 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-N0 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.

Molecules 2018, 23, 2192

24 of 37

Molecules 2018, 23, x

24 of 37

Figure 8. The quantum chemical mechanism of Homo sapiens AChE inhibition by atrazine. The Figure 8. The quantum chemical mechanism of Homo sapiens AChE inhibition by atrazine. The extracted extracted geometry of TS1 (a); TS2 (b); TS3 (c); TS4 (d). geometry of TS1 (a); TS2 (b); TS3 (c); TS4 (d).

The atrazine-based atrazine-based hAChE hAChE inhibition inhibition further further flows flows towards towards the the deprotonation deprotonation of of m-N′ m-N0 atom atom The (Figure 8b and Scheme 1) where hTyr341 phenoxy anion serves as a HB-acceptor. The deprotonation (Figure 8b and Scheme 1) where hTyr341 phenoxy anion serves as a HB-acceptor. The deprotonation concerted via viathe theTS2, TS2,ininwhich which the distance between electronegative atoms amounts Å. concerted the distance between thethe electronegative atoms amounts 3.2113.211 Å. The The similar transition states (this time the TS1s) were also observed for propazine and simazine similar transition states (this time the TS1s) were also observed for propazine and simazine (Figure (Figure S43b and Figure respectively), for peculiar pesticides this wasthe therate-limited rate-limited reaction reaction S43b and Figure S44b, S44b, respectively), and and for peculiar pesticides this was (Figure S45a,b). The calculated calculated activation activation energy energy barriers barriers for for the the generation generation of of propazine-hThr83 propazine-hThr83 and and (Figure S45a,b). The simazine-hThr83 TS1 complexes were 8.5 and 8.9 kcal/mol, respectively (Figure S45a and Figure S45b, simazine-hThr83 TS1 complexes were 8.5 and 8.9 kcal/mol, respectively (Figure S45a and Figure S45b, respectively). To atrazine-hTyr341 TS2TS2 geometry (Figure 8b and 1), atrazine is slightly respectively). Toadopt adoptthe the atrazine-hTyr341 geometry (Figure 8bScheme and Scheme 1), atrazine is in-plane rotated towards hTyr341, whereas propazine and simazine suffered more severe longitudinal slightly in-plane rotated towards hTyr341, whereas propazine and simazine suffered more severe movement (Figure S43b and Figure S44b, respectively). Considering that the energy difference longitudinal movement (Figure S43b and Figure S44b, respectively). Considering that thebetween energy the atrazine-hTyr341 reactant state and TS2 is onlystate 1.6 kcal/mol (Figure 9),kcal/mol an extra stabilization difference between the atrazine-hTyr341 reactant and TS2 is only 1.6 (Figure 9), ancomes extra with the formation of atrazine-hTyr341 TS2. The length of atrazine-hTyr124 HB characterizes it stabilization comes with the formation of atrazine-hTyr341 TS2. The length of atrazine-hTyr124 as HBa moderate one. The second level of QM optimization also confirmed the meaningfulness of TS2. As characterizes it as a moderate one. The second level of QM optimization also confirmed the previous, the meaningfulness of TS2 isthe confirmed since theof56% described thethe proton meaningfulness of TS2. As previous, meaningfulness TS2ofisIRC confirmed since 56% gliding of IRC from TS3 to hTyr341, whereas the 44% or IRC illustrated the proton transfer form TS2 for atrazine m-N0 described the proton gliding from TS3 to hTyr341, whereas the 44% or IRC illustrated the proton atom. Similar noticed foratom. propazine-hTyr341 and simazine-hTyr341 TS1 (Figure S43b and and transfer form trends TS2 forwere atrazine m-N′ Similar trends were noticed for propazine-hTyr341 Figure S44b, respectively). Nevertheless, duringS44b, the rotation of atrazine (the movement propazine simazine-hTyr341 TS1 (Figure S43b and Figure respectively). Nevertheless, duringofthe rotation or simazine), themovement TS1 (IS1) isof sustained (Figure S43b and the Figure of atrazine (the propazine or simazine), TS1S44b, (IS1) respectively). is sustained (Figure S43b and

Figure S44b, respectively).

Molecules 2018, 23, 2192 Molecules 2018, 23, x

25 of 37 25 of 37

Scheme 1. The mechanism of Homo sapiens AChE inhibition by atrazine, propazine, and simazine,

Scheme 1. The mechanism verified by QM studies. of Homo sapiens AChE inhibition by atrazine, propazine, and simazine, verified by QM studies. the third of the reactionpathway, pathway,the the electrophilic electrophilic nature ofof thethe re-established secondary In theInthird stepstep of the reaction nature re-established secondary m-N′ amine is expressed by the release of proton towards the hTyr124 hydroxyl group anion. The 0 m-N amine is expressed by the release of proton towards the hTyr124 hydroxyl group anion. The proton proton transfer takes place via the TS3, in which the distance between the electronegative atoms transfer takes place via the TS3, in which the distance between the electronegative atoms amounts amounts 3.522 Å (Figure 8c and Scheme 1). For the purpose of TS3 generation, atrazine is forced to 3.522 move Å (Figure 8c and Scheme of TS3 generation, is forced to move longitudinally in order1). to For formthe thepurpose tetrahedral structure with theatrazine HB-involved residues. longitudinally in order to form the tetrahedral structure with the HB-involved residues. Therefore, Therefore, this step is slightly energetically more expensive, with the free energy difference of 4.8 this step iskcal/mol slightly(Figure energetically more expensive, the freeTS3 energy difference 4.8 kcal/mol (Figure 9). 9). The characteristic of thewith tetrahedral structure is that of it holds the negatively charged secondary nitrogenTS3 atom. The distinct state is supported by the aromatic of them-N0 The characteristic of them-N′ tetrahedral structure is that it holds the negatively chargednature secondary

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-N0 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

Molecules 2018, 23, x

26 of 37

atrazine’s 1,3,5-triazine ring, which can accept the m-N′ amine negative charge by the intramolecular Molecules 2018, 23, 2192 26 of 37 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 glidingtransfer from TS3 to hTyr124, whereasm-N′ the 43% or On the the IRCother illustrated the propazine-hTyr124 proton transfer from TS3 to proton from TS3 to atrazine’s atom. hand, the and the 0 atom. On the other hand, the propazine-hTyr124 and the simazine-hTyr124 interactions atrazine’s m-N simazine-hTyr124 interactions are characterized by the IS2 (Figure S43c and Figure S44c, are characterized by thethe IS2QM (Figure S43c and Figure S44c, respectively), in which thepropazine QM optimized respectively), in which optimized distances between the hTyr124 and the and distances between the hTyr124 and the propazine and simazine were 4.882 and 3.862 Å, respectively. simazine were 4.882 and 3.862 Å, respectively. Moreover,the theformation formationofof atrazine-hAChE tetrahedral structure appears the Moreover, thethe atrazine-hAChE tetrahedral TS3TS3 structure appears to beto thebeprepre-condition for hHis447 be jointly arranged the atrazine (Figure 8d Scheme and Scheme In last this condition for hHis447 to beto jointly arranged with with the atrazine (Figure 8d and 1). In1). this last step of the reaction pathway, the atrazine’s m-N” amine suffers the nucleophilic attack from step of the reaction pathway, the atrazine’s m-N″ amine suffers the nucleophilic attack from hHis447 hHis447 imidazole ring,the in which proton transferout is carried out from the pesticide the imidazole ring, in which protonthe transfer is carried from the pesticide towards the towards amino acid amino acid via the TS4 (the HB distance was 2.699 Å). Consequently, there is a localization of the via the TS4 (the HB distance was 2.699 Å). Consequently, there is a localization of the negative charge negative around m-N” amine asonce well,again which can be delocalized once again by efficiently delocalized around thecharge m-N″ amine asthe well, which can be efficiently the atrazine′s 1,3,50 by the atrazine s 1,3,5-triazine ring. The similar pathway is observed for the propazine-hHis447 triazine ring. The similar pathway is observed for the propazine-hHis447 and simazine-hTyr447 and simazine-hTyr447 interaction, with the transition stateHB labelled as TS2 (the2.583, HB distances interaction, with the transition state labelled as TS2 (the distances were and 2.540were Å, 2.583, and 2.540 Å, respectively). In comparison with the TS4 (TS2) optimized MD snapshot, the respectively). In comparison with the TS4 (TS2) optimized MD snapshot, the orientation of hHis447 orientation hHis447 needs be adjusted very slightly, there transfer is a proton transfer needs to beof adjusted very toslightly, following which following there is which a proton from the from the atrazine/propazine/simazine m-N” amine to the residue accompanied with the spontaneous atrazine/propazine/simazine m-N″ amine to the residue accompanied with the spontaneous break of break of the scissile m-N”-H bond 8d, Figure S43dS44d, and Figure S44d, respectively). As a the scissile m-N″-H bond (Figure 8d, (Figure Figure S43d and Figure respectively). As a result, a short result, a short and low-barrier hydrogen bond (LBHB) is formed between the hHis447 and the 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), atrazine’s/propazine’s/simazine’s m-N″ amine (Figure 8d, Figure S43d and Figure S44d, and the willing proton transferproton is thus observed. thisEnergetically, is the most favourable reaction respectively), and the willing transfer is Energetically, thus observed. this is the most since the free energy differences between the reactant and transition states were 0.70, 0.68, and favourable reaction since the free energy differences between the reactant and transition states were 0.62 0.68, kcal/mol, respectively. Alongside with the TS1, thethe TS4 (TS2) by (TS2) far the most 0.70, and 0.62 kcal/mol, respectively. Alongside with TS1, theisTS4 is by farimportant the most transition state as it interrupts the hHis447 to serve as a base for hAChE catalytic triad. By By all important transition state as it interrupts the hHis447 to serve as a base for hAChE catalytic triad. means, the meaningfulness of TS4 is confirmed since the 51% of IRC described the proton gliding all means, the meaningfulness of TS4 is confirmed since the 51% of IRC described the proton gliding fromTS3 TS3for forhHis447, hHis447,whereas whereasthe the49% 49%or orIRC IRC illustrated illustratedthe theproton protontransfer transferform formTS4 TS4for foratrazine atrazine from 0 atom (Figure 8d and Scheme 1). Similar trends were also noticed for propazine-hTyr341 and m-N m-N′ atom (Figure 8d and Scheme 1). Similar trends were also noticed for propazine-hTyr341 and simazine-hTyr341 TS2 S43d andand Figure S44d, S44d, respectively). Nevertheless, during the formation simazine-hTyr341 TS2(Figure (Figure S43d Figure respectively). Nevertheless, during the of the atrazine-hHis447 TS4 (propazine-hHis447 and simazine-hTyr447 TS2), the TS1 (IS1), TS2 (TS1), formation of the atrazine-hHis447 TS4 (propazine-hHis447 and simazine-hTyr447 TS2), the TS1 (IS1), and(TS1), TS3 (IS2) sustained. TS2 andwere TS3 (IS2) were sustained.

Figure forfor Homo sapiens AChE inhibition by atrazine determined by B3LYP (6Figure9.9.Free Freeenergy energyprofile profile Homo sapiens AChE inhibition by atrazine determined by B3LYP 31G*) QMQM simulations. (6-31G*) simulations.

The quantum mechanical studies do not provide the basis for the hSer203 covalent modification 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 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 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 particular findings since the metabolites from atrazine, simazine and propazines are not able to

Molecules 2018, 23, 2192

27 of 37

Molecules 2018, 23, x

27 of 37

establish establish the the interactions interactions with with acetylcholinesterase acetylcholinesterase [83]. [83]. This This statement statement was was further further confirmed confirmed in in the the concentration dependent kinetic studies. concentration dependent kinetic studies. 2.8. The Concentration-Dependent Kinetics of Human Acetylcholinesterase Inhibition by 2.8. The Concentration-Dependent Kinetics of Human Acetylcholinesterase Inhibition by 2-Chloro-1,3,5-triazine-based Pesticides 2-Chloro-1,3,5-triazine-based Pesticides There are numerous evidences that organophosphates are irreversible AChE inhibitors whereas 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 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 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 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 related compounds individually influence the hAChE catalytic activity is limited. Therefore, the 22-chloro-1,3,5-triazine-based pesticides are herein evaluated as hAChE inhibitors, in the concentrations 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). 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 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 sapiens AChE activity in a concentration-dependent manner (Figure 10). The inhibition curves implies that the aging reaction, in which the hAChE releases itself from covalently bound pesticide, analysis implies that the aging reaction, in which the hAChE releases itself from covalently bound occurs in the presence of atrazine, propazine, and simazine at a high rate, and that the regarded pesticide, occurs in the presence of atrazine, propazine, and simazine at a high rate, and that the pesticides most likely do no inhibit the hAChE covalently [85]. Therefore, the protocol by which the regarded pesticides most likely do no inhibit the hAChE covalently [85]. Therefore, the protocol by Homo sapiens LD50 values were calculated, as well as all the performed molecular modelling studies of which the Homo sapiens LD50 values were calculated, as well as all the performed molecular modelling 2-chloro-1,3,5-triazine-based pesticides, were experimentally validated. studies of 2-chloro-1,3,5-triazine-based pesticides, were experimentally validated.

Figure 10. 10. Human Human AChE AChE inhibition inhibition by by 2-chloro-1,3,5-triazine-based 2-chloro-1,3,5-triazine-based pesticides. pesticides. Figure

3. Materials and Methods 3. Materials and Methods 3.1. Materials Materials 3.1. Acetylcholine chloride chloride (CAS: (CAS: 60-31-1), 60-31-1), acetylcholinesterase acetylcholinesterase from from human human erythrocytes erythrocytes (CAS: (CAS: Acetylcholine 0 -dithiobis(2-nitrobenzoic acid)) 9000-81-1), and and DTNB DTNB (5,5 (5,5′-dithiobis(2-nitrobenzoic acid)) (CAS: (CAS: 69-78-3), 69-78-3), atrazine atrazine (CAS: (CAS: 1912-24-9), 1912-24-9), 9000-81-1), propazine (CAS: 139-40-2), simazine (CAS: 122-34-9), and PBS (MDL: MFCD00131855) were propazine (CAS: 139-40-2), simazine (CAS: 122-34-9), and PBS (MDL: MFCD00131855) were purchased purchased from Sigma Aldrich (St. Louis, Mo., USA). from Sigma Aldrich (St. Louis, Mo., USA). 3.2. The The Generation Generation of of Pesticides Pesticides Structures Structures 3.2. The crystal crystalstructures structuresofofsimazine simazine monocrotophos [10], dimethoate acetamiprid The [9],[9], monocrotophos [10], dimethoate [11],[11], and and acetamiprid [12] [12] were retrieved Cambridge Crystallographic Data Centre (CCDC)(CSD (CSDIDs: IDs: KUYXIM, were retrieved fromfrom TheThe Cambridge Crystallographic Data Centre (CCDC) ULEJIF,IPCPYB, IPCPYB, HANBAA). HANBAA). The The remaining remaining pesticides pesticides structures structures were were generated generated by by drawing drawing using using ULEJIF, the MacroModel MacroModel version version 88 software software(Schrodinger, (Schrodinger,Cambridge, Cambridge,MA, MA,USA) USA)[75]. [75]. the 3.3. 3.3. The The Conformational Conformational Analysis Analysis The The conformational conformational analyses analyses were were performed performed by by using using the the MacroModel MacroModel version version 88 software software [75] [75] via approach [76]. [76]. The solute was described via the theMonte MonteCarlo/multiple Carlo/multiple minimum minimum (MC/MM) (MC/MM) approach described using using five different potentials (MM3, AMBER94, MMFF, MMFFs, OPLSAA) [77,85–88]. The effect of solvent was incorporated into the MC/MM calculations using the generalized Born/surface area (GB/SA)

Molecules 2018, 23, 2192

28 of 37

five different potentials (MM3, AMBER94, MMFF, MMFFs, OPLSAA) [77,85–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

Molecules 2018, 23, 2192

29 of 37

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,

Molecules 2018, 23, 2192

30 of 37

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): xA = f rmsd, ≤ a + 0.5( f rmsd ≤ b − f rmsd ≤ a) (7) In particular: xA = CAA or FFDAA, in the case of LB alignment accuracy, whereas the xA = DA, in the case of SB accuracy. The frmsd ≤ a and frmsd ≤ b 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.

Molecules 2018, 23, 2192

31 of 37

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): ∆Gbind = ∆Eele + ∆EvdW + ∆Gsolv + T∆S,

(8)

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.

Molecules 2018, 23, 2192

32 of 37

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ˇca, 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.

Molecules 2018, 23, 2192

33 of 37

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

References 1. 2.

3.

4. 5. 6. 7.

8.

9. 10.

11.

12. 13. 14. 15.

16. 17. 18. 19. 20. 21.

Gilden, R.C.; Huffling, K.; Sattler, B. Pesticides and health risks. J. Obstet. Gynecol. Neonatal Nurs. 2010. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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). Bolton, T.B.; Lim, S.P. Action of acetylcholine on smooth muscle. Z. Kardiol. 1991, 80 (Suppl. 7), 73–77. [PubMed] Jones, B.E. From waking to sleeping: Neuronal and chemical substrates. Trends Pharmacol. Sci. 2005, 26, 578–586. [CrossRef] [PubMed] Chapalamadugu, S.; Chaudhry, G.R. Microbiological and biotechnological aspects of metabolism of carbamates and organophosphates. Crit. Rev. Biotechnol. 1992. [CrossRef] [PubMed] Đorđevi´c, J.; Vladisavljevi´c, G.T.; Trti´c-Petrovi´c, 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. [CrossRef] [PubMed] Đorđevi´c, J.S.; Vladisavljevi´c, G.T.; Trti´c-Petrovi´c, 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. [CrossRef] 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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. [CrossRef] 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. [CrossRef] Boudina, A.; Emmelin, C.; Baaliouamer, A.; Grenier-Loustalot, M.F.T.; Chovelon, J.M. Photochemical behaviour of carbendazim in aqueous solution. Chemosphere 2003. [CrossRef] 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. 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] Park, C. A Dictionary of Environment and Conservation, 1st ed.; Oxford University Press: Oxford, UK, 2007. 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. [CrossRef] [PubMed] USEPA. Reregistration Eligibility Decision (RED), Linuron. 1995. Available online: https://archive.epa.gov/ pesticides/reregistration/web/pdf/0047.pdf (accessed on 15 June 2018). 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. 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).

Molecules 2018, 23, 2192

22. 23. 24. 25. 26.

27. 28. 29.

30. 31.

32. 33. 34.

35.

36.

37. 38. 39.

40. 41.

42. 43.

34 of 37

Extension Toxicology Network. Pesticide Information Profiles, Oregon State University, Diuron. 1996. Available online: http://www.extoxnet.orst.edu (accessed on 19 March 2016). Worthing, C.R.; Walker, S.B. The Pesticide Manual: A World Compendium, 8th ed.; The British Crop. Protection Council: Thornton Heath, UK, 1987. 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). Hayes, W.J.J.; Laws, E.R.J. Handbook of Pesticide Toxicology; Academic Press: New York, NY, USA, 1991. 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). Occupational Health Services, Inc. MSDS for Dimethoate. OHS Inc., Secaucus, NJ, USA. 1991. Available online: https://hr.umich.edu (accessed on 17 March 2016). 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). 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). Kidd, H.; James, D. Agrochemicals Handbook, 3rd ed.; Royal Society of Chemistry: Cambridge, UK, 1991. 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). Reagan-Shaw, S.; Nihal, M.; Ahmad, N. Dose translation from animal to human studies revisited. FASEB J. 2008, 22, 659–661. [CrossRef] [PubMed] 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. 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. [CrossRef] [PubMed] 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. [CrossRef] 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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. [CrossRef] 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. [PubMed] 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed]

Molecules 2018, 23, 2192

44.

45.

46.

47.

48. 49.

50. 51. 52.

53. 54. 55. 56.

57. 58. 59. 60. 61. 62. 63. 64. 65.

35 of 37

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). 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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. [CrossRef] 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] Kennedy, G.L. Acute toxicity studies with oxamyl. Fundam. Appl. Toxicol. 1986, 6, 423–429. [CrossRef] Turusov, V.; Rakitsky, V.; Tomatis, L. Dichlorodiphenyltrichloroethane (DDT): Ubiquity, persistence, and risks. Environ. Health Perspect. 2002, 110, 125–128. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] Phosmet. Available online: http://pmep.cce.cornell.edu/profiles/extoxnet/metiram-propoxur/phosmetext.html (accessed on 7 April 2018). Sanderson, D.M. Treatment of poisoning by anticholinesterase insecticides in the rat. J. Pharm. Pharmacol. 1961, 13, 435–442. [CrossRef] [PubMed] 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). 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). Extoxnet Pip Terbufos. Available online: http://extoxnet.orst.edu/pips/terbufos.htm (accessed on 4 July 2018). Dichlorvos NIOSH Publications and Products. Available online: https://www.cdc.gov/niosh/idlh/62737. html (accessed on 2 July 2018). 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). Diazinon. Available online: http://www.herbiguide.com.au/Descriptions/hg_DIAZINON.htm (accessed on 2 July 2018). 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. Pearce, E.M. Kirk-OTHMER encyclopedia of Chemical Technology, 3rd ed.; Wiley-Interscience: New York, NY, USA, 1978. Extoxnet Pip Oxamyl. Available online: http://extoxnet.orst.edu/pips/oxamyl.htm (accessed on 4 July 2018). “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). 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).

Molecules 2018, 23, 2192

66.

67. 68. 69.

70.

71. 72. 73. 74.

75.

76. 77. 78.

79.

80. 81.

82.

83.

84. 85.

36 of 37

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). Methyl Parathion. Available online: http://pmep.cce.cornell.edu/profiles/extoxnet/haloxyfopmethylparathion/methyl-parathion-ext.html#2 (accessed on 2 July 2018). Dicamba (Banvel) Herbicide Profile 10/83. Available online: http://pmep.cce.cornell.edu/profiles/herbgrowthreg/dalapon-ethephon/dicamba/herb-prof-dicamba.html (accessed on 4 July 2018). 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). 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. [CrossRef] [PubMed] 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). 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). ˇ Colovi´ c, M.B.; Krsti´c, D.Z.; Lazarevi´c-Pasti, T.D.; Bondzi´c, AM.; Vasi´c, V.M. Acetylcholinesterase Inhibitors: Pharmacology and toxicology. Curr. Neuropharmacol. 2013. [CrossRef] 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. [CrossRef] [PubMed] 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. [CrossRef] Chang, G.; Guida, W.C.; Still, W.C. An internal-coordinate Monte Carlo method for searching conformational space. J. Am. Chem. Soc. 1989. [CrossRef] 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. [CrossRef] Mladenovi´c, 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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. [CrossRef] Rosenfeld, C.A.; Sultatos, L.G. Concentration-dependent kinetics of acetylcholinesterase inhibition by the organophosphate paraoxon. Toxicol. Sci. 2006. [CrossRef] [PubMed]

Molecules 2018, 23, 2192

86. 87. 88.

89. 90.

91. 92.

93.

94.

95.

96.

97. 98. 99.

100.

101. 102. 103.

37 of 37

Halgren, T.A. Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94. J. Comput. Chem. 1996. [CrossRef] 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. [CrossRef] [PubMed] 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. [CrossRef] 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. [CrossRef] Arsic, B.; Aguilar, J.A.; Bryce, R.A.; Barber, J. Conformational study of tylosin A in water and full assignments of 1 H and 13 C spectra of tylosin A in D2 O and tylosin B in CDCl3 . Magn. Reson. Chem. 2017. [CrossRef] [PubMed] Clark Still, W.; Tempczyk, A.; Hawley, R.C.; Hendrickson, T. Semianalytical treatment of solvation for molecular mechanics and dynamics. J. Am. Chem. Soc. 1990. [CrossRef] 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] 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. 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. [CrossRef] [PubMed] 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. [CrossRef] [PubMed] Resource for Biocomputing, Visualization, and Informatics (RBVI). Available online: http://rbvi.ucsf.edu (accessed on 6 August 2018). 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. 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. [CrossRef] [PubMed] Dennington, R.; Keith, T.; Millam, J. GaussView 2009. Available online: http://gaussian.com/glossary/g09/ (accessed on 13 March 2018). 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. Ellman, G.L.; Courtney, K.D.; Andres, V.; Featherstone, R.M. A new and rapid colorimetric determination of acetylcholinesterase activity. Biochem. Pharmacol. 1961. [CrossRef]

Sample Availability: Samples of the compounds are not available from the authors. © 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).