Binding organophosphate pesticides to acetylcholinesterase: risk assessment using the Monte Carlo method

Abstract

Binding to acetylcholinesterase (AChE k1) may cause toxic effects in humans. Organophosphates simultaneously are both dangerous and useful substances. Dangerous since they are employed in chemical warfare and useful when they are applied as pesticides. Here, we suggest the models for organophosphates binding to AChE k1 developed via representing the molecular structure by a simplified molecular input-line entry system using so-called optimal descriptors calculated with the Monte Carlo technique using the Correlation and Logic (CORAL) free software available on the Internet (http://www.insilico.eu/coral). Quantitative structure-activity relationships (QSARs) serve to develop predictive models for organophosphates. The predictive potential of these models is quite good: the determination coefficient for the validation set ranged from 0.87 to 0.90. These models were built up according to the principle ‘QSARs is a random event’, that is, predictive potential of an approach should be checked up with several splits of available data into the training and test sets. The special scheme of mechanistic interpretation definition is represented. The mechanistic interpretation is based on probabilities of molecular features to be in the sub-group of promoters of increase for endpoint or in sub-group of promoters of it is decrease.

Keywords:

• Organophosphate
• pesticide toxicity
• QSAR
• risk assessment
• SMILES
• Monte Carlo method

Introduction

Pesticides are widely used in agriculture for plant protection and to increase the yields and quality of agricultural products, but also in domestic applications (Naik, Singh, and Singh 2009). Organophosphate compounds (OPC) are the widely used choice as insecticides and are widely used across the world. For this reason, these substances are frequently encountered in surface water and groundwater due to their agricultural use. OPC are also used for ornamentals, in floriculture, and for residential treatments. Unfortunately, despite these advantages, these pesticides are highly toxic and also bio accumulate (Torres-Palma and Serna-Galvis 2018; Bora et al. 2020). Exposure to OPCs can affect the immune response and carbohydrate metabolism, and cause cardiovascular toxicity (Marimuthu et al. 2019).

Residues of these pesticides usually build up on the surface of the soil and ultimately find their way into food products. They may have adverse impacts on various groups of soil microbes, wild animals, and also humans, through environmental contamination, or residues in food and feedstuff (Torres-Palma and Serna-Galvis 2018; Bora et al. 2020; Toropov et al. 2017, 2018). Experimental data on these chemicals can be used to improve the design of new pesticides, moving toward safer, greener substances, considering their adverse effects on humans and ecosystems (Mudasir, Wibowo, and Pranowo 2013; Wibowo, Mudasir, and Pranowo 2017). Furthermore, OPC are chemical warfare tools (Mendonca et al. 2020).

Quantitative structure-activity relationships (QSARs) serve to develop predictive models for the toxicity of compounds in general and OPC in particular (Niraj, Saini, and Kumar 2015; Marimuthu et al. 2019; Bora et al. 2020; Maxwell and Brecht 1992; Ruark et al. 2011; Zhao and Yu 2013; Morrill et al. 2015; Lee and Barron 2016). The Correlation and Logic (CORAL) software (http://www.insilico.eu/coral) has been used up as a tool to develop QSAR models for different endpoints (Toropov et al. 2017, 2018; Veselinović et al. 2015; Kumar and Chauhan 2017; Kumar, Kumar, and Sindhu 2019; Islam and Pillay 2016). Here, we studied the ability of this software to build up predictive models for toxicity OPC binding to acetylcholinesterase (AChE k1).

Method

Data

A database of 278 molecular structures together with numerical data on toxicity represented by rate constants for human binding to AChE k1 were taken from Ruark et al. (2013). These substances have been randomly distributed into an active training set (≈25%) used to construct the model, a passive training set (≈25%) independent of those used for building the model, a calibration set (≈25%) to identify situations under which the model is over-training, and a validation set (≈25%) to provide the final assessment of the statistical quality of the model (Toropov, Toropova, and Benfenati 2020). The molecular structure was represented with simplified molecular input-line entry systems (SMILES; Weininger 1988).

Optimal SMILES-based descriptors

The optimal SMILES-based descriptor of the correlation weights (DCW) applied here is calculated as follows: (1) $${\mathit{DCW}}_{\mathit{SMILES}}\left(T^{\mathrm{*}},N^{\mathrm{*}}\right)=CW\left(\mathit{APP}\right)+\sum _{\mathrm{ }}^{\mathrm{ }}CW(S_k)+\sum _{\mathrm{ }}^{\mathrm{ }}CW({SS}_k)$$(1)

Correlation weights (CW) are calculated with the Monte Carlo method using the CORAL software (http://www.insilico.eu/coral). T is an integer to define rare fragments and active ones. N is the number of epochs of the Monte Carlo optimization. Atom-pairs proportion (APP) is atom pair’s proportions self-organized vector.

In Table 1, examples of APP are given. S_k is SMILES atoms, that is, one symbol or a group of symbols that cannot be examined separately (‘Br’, ‘%11’, ‘@@’). SS_k is fragments of SMILES containing two SMILES atoms. CW(APP), CW(S_k), and CW(SS_k) are correlation weights for these SMILES attributes.

Table 1. An example of building a self-organized vector of APP for the following SMILES: CCC(Br)(CC)C(=O)NC(N)=O.

Atom 1	Atom 2	n1	n2	Code applied to build up the model
Br	N	1	2	(Br.N)..1.2.
Br	O	1	2	(Br.O)..1.2.
Br	Double bond	1	2	(Br.=)..1.2.
N	O	2	2	(N..O)..2.2.
N	Double bond	2	2	(N..=)..2.2.
O	Double bond	2	2	(O..=)..2.2.

T* is the threshold, that is, a value to discriminate molecular features extracted from SMILES into two categories, viz. rare (frequency in the active training set is less than T*) and active (frequency in the active training set is equal or better than T*). N* is the number of epochs of the Monte Carlo optimization.

Monte Carlo optimization

QSAR-models were calculated with the Monte Carlo optimization based on the target functions: (2) $${TF}_{}=r_{AT}+r_{PT}-\left|r_{AT}-r_{PT}\right|\mathrm{*}0.1$$(2)

$r_{AT}$ and $r_{PT}$ are correlation coefficients between observed and predicted endpoint for the active training and passive training sets, respectively. The maximum given by Equation 2(2) $${TF}_{}=r_{AT}+r_{PT}-\left|r_{AT}-r_{PT}\right|\mathrm{*}0.1$$(2) is the target of the Monte Carlo calculations. Having the corresponding numerical data in the CWs, the endpoint is calculated as follows: (3) $$\mathit{human}\mathrm{ }\mathit{AChE}\mathrm{ }{K}_1=C_0+C_1{\times \mathit{DCW}}_{\mathit{SMILES}}\left(T^{\mathrm{*}},N^{\mathrm{*}}\right)$$(3) the C_O and C₁ are regression coefficients.

Results and discussion

Domain of applicability

The domain of applicability for these compounds is defined according to so-called SMILES defects (Toropova et al. 2020). The number of potential outliers in the validation set is 6, 7, and 5, respectively, for splits 1, 2, and 3. The SMILES defect is not a rigid categorization of compounds but only indicates that the representability of fragments of the corresponding SMILES in the active training set (Toropova et al. 2020) is not enough. Figure 1 gives examples of the statistical outliers. The outliers are not all similar for different splits though some are the same (e.g. #1 for splits 1 and 2; #198 for splits 2 and 3).

Figure 1. Statistical outliers for the validation set in the cases split 1, 2, and 3.

Mechanistic interpretation

In Table 2, promoters (Toropova et al. 2020) for increases and decreases for rate constants of human AChE k1 are listed. There are some promoters of an increase (double bonds indicated by ‘=’) since in all three runs of the Monte Carlo optimization this SMILES attribute obtains positive CWs (Table 2) and some promoters of a decrease (brackets in SMILES) since this SMILES attribute obtains negative CWs in all three runs of the optimization. However, there are some more complex and more intriguing promoters, for example, oxygen and phosphorus in proportion 2:1 [SMILES attribute “(O.P).2.1.”] is a promoter of greater toxicity, while the fragment PO [SMILES attribute “P…O…….”] is a promoter of lower toxicity. These data can be useful from the heuristic point of view. It is to be noted that (O.P).2.1. means the presence oxygen and phosphoros atoms in ration 2:1 (maybe not connected directly); whereas P…O……. means connecting of these atoms by covalent bond.

Table 2. Mechanistic interpretation – promoters of increase and decrease for the toxicity of the pesticides.

SMILES-attribute SAk	CWs* run 1	CWs run 2	CWs run 3	N1	N2	N3	Defect dk
Promoters of an increase
=………..	0.96519	0.42842	0.42862	71	69	69	0.0000
C…(…….	0.06438	0.05487	0.08783	71	69	69	0.0000
P………..	5.25534	6.58765	6.14333	71	69	69	0.0000
C…C…….	0.33947	0.03825	0.14633	67	67	66	0.0001
O…C…….	0.71139	0.51058	0.90372	64	62	67	0.0005
(P..=)..1.1.	1.11519	0.19485	0.41416	62	56	47	0.0018
=…(…….	0.69470	0.75006	0.41743	37	35	42	0.0011
(O..P)..2.1.	0.36520	0.50449	0.87799	31	24	24	0.0016
(O..=)..2.1.	1.25885	1.67468	1.13817	28	20	15	0.0041
(N..P)..1.1.	1.45260	0.89589	0.28152	27	26	27	0.0002
(O..P)..3.1.	0.41247	1.16865	0.81350	26	27	31	0.0015
(O..=)..3.1.	1.87320	1.78687	0.30757	25	27	27	0.0008
[S+]……..	4.70518	4.98571	4.01479	16	8	7	0.0054
N…C…….	0.67927	0.53490	0.32812	13	9	11	0.0010
(S..P)..2.1.	0.98399	0.75977	1.06371	12	12	20	0.0038
Promoters of a decrease
(………..	–0.02859	–0.32915	–0.07873	71	69	69	0.0000
P…(…….	–3.36380	–4.03716	–2.39460	71	69	69	0.0000
S………..	–0.60222	–0.69373	–0.16600	62	59	62	0.0002
S…C…….	–0.18549	–0.56536	–0.70620	62	59	62	0.0002
1………..	–0.66790	–0.66383	–0.93290	56	49	41	0.0020
(S..=)..1.1.	–1.08107	–0.23634	–1.12585	47	43	37	0.0015
P…=…….	–0.10617	–0.41484	–0.20683	35	36	29	0.0011
C…1…….	–0.43960	–0.37729	–0.12193	23	16	12	0.0043
N………..	–0.42554	–0.45419	–0.12301	20	18	20	0.0002
1…(…….	–0.40147	–0.81489	–0.33245	14	15	12	0.0009
C…2…….	–0.61325	–0.81121	–0.17312	13	9	9	0.0024
P…O…….	–0.38286	–0.61044	–0.44162	13	20	19	0.0029
(P..=)..1.2.	–1.27212	–1.31094	–0.98613	8	9	19	0.0060
2…(…….	–0.57855	–0.26338	–0.15654	7	13	4	0.0037
N…(…….	–0.03706	–0.82781	–1.04743	6	7	5	0.0011

*CWs, correlation weights; N1, N2, and N3 are frequencies of a SMILES attribute in the active training set, passive training set; and calibration set, respectively. The d_k is the statistical defect of SMILES attribute, calculated via probability in active training set P(SA_k), probability in calibration set P′(SA_k); frequencies in active training set N(SA_k), in calibration set N′(SA_k), as the following:

$d_k=\frac{\left|{P(SA}_k)-{P\prime (SA}_k)\right|}{N\left({SA}_k\right)+N\prime ({SA}_k)}$

Models

Computational experiments with three random splits give the following models:

Split 1

(4) $$\mathit{human} \mathit{AChE} K_1=-0.9042(\pm \mathrm{ }0.0411)\mathrm{ }+\mathrm{ }0.6609(\pm \mathrm{ }0.0043)\text{ * DCW}(1,15)$$(4)

Split 2

(5) $$\mathit{human} \mathit{AChE} K_1=1.0303(\pm \mathrm{ }0.0275)\mathrm{ }+\mathrm{ }1.0200(\pm \mathrm{ }0.0063)\text{ * DCW}(1,15)$$(5)

Split 3

$\mathit{human}\mathrm{ }\mathit{AChE}\mathrm{ }{K}_1=0.9458(\pm \mathrm{ }0.0239)+0.8805(\pm \mathrm{ }0.0050)\mathrm{*}\mathrm{DCW}(1,15)$(6)

Table 3 contains the statistical characteristics of models for three random splits. There is dispersion of the statistical characteristics, but not too much. In other words, the models suggested are robust and can be used for practical applications.

Table 3. Statistical characteristics of models for splits 1–3.

Split	Set	N*	R^2	Q^2	RMSE	F
1	Active training	71	0.7670	0.7557	0.868	227
	Passive training	69	0.8756	0.8697	0.786	472
	Calibration	69	0.8913	0.8858	0.577	549
	Validation	69	0.8700		0.568
2	Active training	71	0.7991	0.7891	0.794	275
	Passive training	69	0.8112	0.8013	0.910	288
	Calibration	69	0.8685	0.8615	0.560	443
	Validation	69	0.9041		0.537
3	Active training	71	0.7908	0.7812	0.856	261
	Passive training	69	0.8193	0.8092	0.781	304
	Calibration	69	0.7449	0.7225	0.890	196
	Validation	69	0.8950		0.590

F: Fischer F-ratio; *N: number of compounds in a set; Q: cross-validated R; R: correlation coefficient; RMSE: root mean squared error.

Comparison with models from the literature

In Table 4, the statistical quality of models suggested here is compared with some other models suggested. Despite the interest in OPCs (Devillers 2000; Mudasir, Wibowo, and Pranowo 2013; Ruark et al. 2013; Bora et al. 2020) it is impossible to find any QSAR studies on the rate constants of human AChE k1. Therefore the only comparison can be between the results here and those from the original article (Ruark et al. 2013). In our case, the original set of compounds was split into four sets: active training, passive training, calibration, and validation, while Ruark et al. (2013) report the results on three subsets. However, it is possible to make a comparison. The results of the active training are quite similar. No results on the passive training are available in the original article, because the authors applied a different approach; however, the values for the passive training used only in our case are higher than those on the training set in the original article. The values for the calibration and validation sets are always better in our case, particularly for the validation set. The validation set values are of particular interest because they demonstrate that our model is not affected by overtraining and quite probably can be applied for substances not used in the modeling phase.

Table 4. Comparison with statistical characteristics of a model described in the literature.

Reference	R^2 active training set	R^2 passive training set	R^2 calibration set	R^2 validation set
Ruark et al. (2013)	0.80	–	0.76	0.66
Here
Split 1	0.767	0.875	0.891	0.870
Split 2	0.800	0.811	0.869	0.904
Split 3	0.791	0.819	0.745	0.895

Besides their good statistical results, our model is simpler, does not require the calculation of molecular descriptors as in Ruark’s article (2013), but only needs the SMILES to represent the molecular structure.

Supplementary materials section gives technical details of the suggested models.

Conclusions

We have developed new models for AChE binding for OPC that are robust and predictive. These models are accompanied by the domain of applicability and a mechanistic interpretation. The statistical quality of the models examined here is better than in models for the same endpoint obtained with the CODESSA software (Ruark et al. 2013). The approach we adopt is quite simple and uses the freely available CORAL software.

Ethical approval

The manuscript is not submitted to any other journal.

Consent to participate

The article did not involve any human participants.

Consent to publish

The submitted work is original.

Author contributions

APT was involved in formal analysis, data collection, and curation, investigation, methodology incl. statistics, software, visualization, writing of original draft, and revision. AAT was involved in conceptualization, data curation, formal analysis, investigation, methodology incl. statistics, software, visualization, writing of original draft, review, and editing. AR was involved in conceptualization, formal analysis, visualization, writing of original draft, review, and editing. EB was involved in conceptualization, formal analysis, visualization, review, and editing.

Acknowledgments

The authors are grateful for the contribution of the project LIFE-CONCERT contract (LIFE17 GIE/IT/000461) and Optitox (OC/EFSA/SCER/2018/01) for support.

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this article.

Data availability statement

All data generated or analyzed during this study are included in this article (and its supplementary information files).

Additional information

Funding

Not applicable.