# Navigating the complexities of docking tools with nicotinic receptors and acetylcholine binding proteins in the realm of neonicotinoids

Zakaria Bouchouireb, Damien Olivier-Jimenez, Titouan Jaunet-Lahary, Steeve H Thany, Jean-Yves Le Questel

<figure><img src="image_1.png" /></figure>

## To cite this version:

Zakaria Bouchouireb, Damien Olivier-Jimenez, Titouan Jaunet-Lahary, Steeve H Thany, Jean-Yves Le Questel. Navigating the complexities of docking tools with nicotinic receptors and acetylcholine binding proteins in the realm of neonicotinoids. Ecotoxicology and Environmental Safety, 2024, 281, pp.116582. 10.1016/j.ecoenv.2024.116582. hal-04824048

HAL Id: hal-04824048

https://hal.science/hal-04824048v1

Submitted on 6 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Navigating the complexities of docking tools with nicotinic receptors and acetylcholine binding proteins in the realm of neonicotinoids

Zakaria Bouchuireb<sup>a, *</sup>, Damien Olivier-Jimenez<sup>b</sup>, Titouan Jaunet-Lahary<sup>a</sup>, Steeve H. Thany<sup>c, d</sup>, Jean-Yves Le Questel<sup>a, *</sup>

<figure><img src="image_4.png" /></figure>

Check for
updates

## ARTICLE INFO

**Keywords:**

Benchmark

Molecular docking

Autodock

Ledock

CDOCKER

AChBPs

nAChRs

Neonicotinoids

## ABSTRACT

Molecular docking, pivotal in predicting small-molecule ligand binding modes, struggles with accurately identifying binding conformations and affinities. This is particularly true for neonicotinoids, insecticides whose impacts on ecosystems require precise molecular interaction modeling. This study scrutinizes the effectiveness of prominent docking software (Ledock, ADFR, Autodock Vina, CDOCKER) in simulating interactions of environmental chemicals, especially neonicotinoid-like molecules with nicotinic acetylcholine receptors (nAChRs) and acetylcholine binding proteins (AChBPs). We aimed to assess the accuracy and reliability of these tools in reproducing crystallographic data, focusing on semi-flexible and flexible docking approaches. Our analysis identified Ledock as the most accurate in semi-flexible docking, while Autodock Vina with Vinardo scoring function proved most reliable. However, no software consistently excelled in both accuracy and reliability. Additionally, our evaluation revealed that none of the tools could establish a clear correlation between docking scores and experimental dissociation constants ($K_d$) for neonicotinoid-like compounds. In contrast, a strong correlation was found with drug-like compounds, bringing to light a bias in considered software towards pharmaceuticals, thus limiting their applicability to environmental chemicals. The comparison between semi-flexible and flexible docking revealed that the increased computational complexity of the latter did not result in enhanced accuracy. In fact, the higher computational cost of flexible docking with its lack of enhanced predictive accuracy, rendered this approach useless for this class of compounds. Conclusively, our findings emphasize the need for continued development of docking methodologies, particularly for environmental chemicals. This study not only illuminates current software capabilities but also underscores the urgency for advancements in computational molecular docking as it is a relevant tool to environmental sciences.

## 1. Introduction

Neonicotinoids, introduced in the latter part of the 20th century, have been heralded as a significant advancement in the domain of insecticides. Their systemic action, allowing them to be absorbed by plants, presented a unique mechanism to tackle a wide spectrum of pests (Jeschke et al., 2011). Their use in agricultural settings was not just due to their efficacy, but also because they emerged as an alternative to older insecticides that were losing the battle against pest resistance (Casida and Durkin, 2013) and appeared more toxic for mammals (Marrs, 2012). However, the benefits of neonicotinoids have been shadowed by

environmental concerns. The extent of their use, combined with their persistence, has led to detections in various ecosystems. Studies have reported residues in soils, waterways, and even in plants that were never directly treated with these insecticides, painting a worrisome picture of environmental contamination (van der Sluijs et al., 2014).

The repercussions of such widespread contamination are not just limited to the targeted pests. There is growing evidence to suggest that neonicotinoids have been adversely affecting non-target organisms. For instance, they have been implicated in the decline of beneficial pollinators. Bees, pivotal to our ecosystems, have been particularly affected, with studies implicating neonicotinoids in their decline. Even sub-lethal

<footnote><sup>a</sup> Nantes Université, CNRS, CEISAM, UMR 6230, Nantes 44000, France</footnote>

<footnote><sup>b</sup> Leiden University Medical Center, Center for Proteomics and Metabolomics, Albinusdreef 2, Leiden 2333ZA, Netherlands</footnote>

<footnote><sup>c</sup> Université d'Orléans, Physiology, Ecology and Environment (P2E) laboratory USC INRAE 1328, Orléans 45067, France</footnote>

<footnote><sup>d</sup> Institut universitaire de France (IUF), 1 rue Descartes 75005 Paris, France</footnote>

<footnote>* Corresponding authors.</footnote>

<footnote>E-mail addresses: zakaria.bouchuireb@univ-nantes.fr (Z. Bouchuireb), jean-yves.le-questel@univ-nantes.fr (J.-Y. Le Questel).</footnote>
doses have been shown to impact bee foraging behavior, reproductive success, and overall colony health (Whitehorn et al., 2012; Henry et al., 2012; Gill et al., 2012; Rundlöf et al., 2015; Pisa et al., 2014; Woodcock et al., 2016; Tsvetkov et al., 2017). Furthermore, concerns about the potential impacts of neonicotinoids on mammalian health are growing, warranting deeper investigation into their effects on human health and safety (Han et al., 2018; Wang et al., 2020; Chen et al., 2020; Zhang et al., 2022).

Neonicotinoids act as competitive modulators, specifically targeting the nicotinic acetylcholine receptors (nAChRs) in insects. With the rising ecological concerns associated with neonicotinoids, computational biology and cheminformatics have stepped in, offering innovative tools to understand and forecast the dynamics of these insecticides across their targets. Among these tools, molecular docking stands out, having been extensively employed as a cheap and scalable method to anticipate molecular interactions, particularly in drug discovery realms (Morris and Lim-Wilby, 2008). Essentially, molecular docking facilitates the prediction of how molecules orient themselves and bind, establishing stable complexes with biological entities, thus shedding light on potential interactions and their subsequent implications. Several docking software options, including Autodock (Morris and Lim-Wilby, 2008), Autodock Vina (Trott and Olson, 2009; Eberhardt et al., 2021), Ledock (Zhao and Caflisch, 2013), UCSF DOCK (Coleman et al., 2013), GOLD (Verdonk et al., 2003) and CDOCKER (Wu et al., 2003; Buckner et al., 2023) are accessible for both academic and industry purposes. These tools rely on a combination of sampling algorithms and scoring functions to generate results. While many of these programs often align well with experimental findings, variations in their algorithms and functions can lead to diverse outputs (Warren et al., 2005; Huang and Zou, 2010; Cross et al., 2009). Nevertheless, these software packages prove invaluable in assessing risks to the environment and human health. For instance, Mansouri et al. employed them to predict the estrogen receptor activity of man-made chemicals, which is crucial for understanding their potential endocrine-disrupting effects (Mansouri et al., 2016). Similarly, Xue et al. leveraged docking techniques to corroborate their findings on the evaluation of binding performance of flavonoids against the same family of receptors (Xue et al., 2022). Other studies, like those by CarterEAU et al. and Alamiddine et al., blended docking with other methodologies like molecular dynamics simulations and patch-clamp recordings to delve into the binding process of neonicotinoids, namely imidacloprid, thiacloprid and sulfoxaflor against insect's nAChRs (CarterEAU et al., 2021; Alamiddine et al., 2019).

However, as we venture into the world of neonicotinoids, challenges emerge. As eluded to by several studies (Schneider and Fechner, 2005; Cross et al., 2009), these computational tools, despite their sophistication, possess inherent biases linked to the target datasets they were optimized for. Tailored predominantly for pharmaceutical research, these tools might be less equipped to accurately handle and predict the interactions of molecules like neonicotinoids. Adding another layer of complexity is their molecular architecture. They often incorporate Pan Assay Interference Compounds (PAINS) substructures, which are moieties generally avoided in pharmaceutical research due to their potential to interfere in assays and give false-positive results (Baeil and Holloway, 2010). Traditional docking software, trained on datasets dominated by pharmaceutical compounds, might be ill-prepared to accurately model and predict the behavior of neonicotinoids, given the presence of these PAINS substructures. Moreover, the benchmarks and datasets that underpin the evaluation and validation of these molecular docking tools often overlook the unique challenges posed by ecotoxicological compounds. As highlighted by Xue et al. (2022), most benchmark datasets are skewed towards pharmaceutical compounds, potentially bypassing the unique characteristics and challenges associated with agricultural and ecotoxicological compounds. This issue is further highlighted by Chaput et al. (2016), where an adjustment of the DUD-E dataset led to a significant reduction in virtual screening success rates specifically for commercial software, as assessed in their study. More recently, similar

biases were observed by Chen et al. (2019), Sieg et al. (2019) in the context of machine learning models developed from classical datasets. Thus, this inherent bias might lead to over-optimistic evaluations of these tools when applied outside their primary domain of pharmaceuticals.

The benchmarks and datasets commonly used for evaluating molecular docking tools often don't cater to the specific challenges presented by compounds relevant to ecotoxicology. As underscored by Cheng et al. (2012) and further discussed by Xue et al. (2022), when it comes to computational toxicology, docking tools face two main challenges. First, they must tackle the issue of low binding affinity. Indeed, for human and environmental risk assessment, it's important to note that many of the chemicals studied exhibit mainly weak interactions with their target sites. Second, they need to accurately rank ligands based on this binding affinity, a task that becomes particularly complex when dealing with chemicals that have closely related structures. Consequently, when analyzing vast datasets, these weakly interacting chemicals may skew the results, behaving as statistical outliers. As such, when evaluating the effectiveness of these programs, it's crucial to clarify their intended application, be it for environmental concerns or human health risk assessments.

Merging the fields of molecular docking and ecotoxicology promises new insights, but remains an emerging area. As the environmental impact of pollutants like neonicotinoids grows, the need for computational methods that can rapidly provide clarity on their interactions, behaviors, and potential countermeasures becomes more pressing. For these computational forecasts to truly guide action, it's essential that the tools undergo thorough evaluation against datasets reflecting the specific ecotoxicological challenges. To this end, we examined the PDB dataset, featuring 54 protein-ligand interaction systems, where proteins are either nAChRs or acetylcholine binding proteins (AChBPs) paired with small molecule ligands from the neonicotinoid family or related compounds. We re-docked the co-crystallized ligands using four notable docking software that are widely used, easy to interface with interpreted languages to write workflows and freely available for academics: (i) Ledock, (ii) Autodock Flexible Receptor (ADFR), (iii) Autodock Vina, and (iv) CDOCKER. Although one could argue that excluding top-rated commercial software is a limitation, this software selection was driven by a commitment to evaluate and discuss tools that are readily available for broader academic and research use, particularly in contexts where access to commercial software might be limited. Our study contrasts the docking accuracy of these prominent software tools in predicting the binding mechanisms of neonicotinoids with nAChRs and AChBP. Moreover, we gauged the precision of these tools in forecasting correct binding affinities and their order, cross-referencing with existing experimental data for select systems. Using a dataset tailored for ecotoxicological concerns, we seek to identify the best docking tool to accurately model the complex interactions of neonicotinoids against their targets, filling the current knowledge void and informing subsequent scientific endeavors and policy decisions.

## 2. Methods

### 2.1. *Systems parsing and initial preparation*

The basis for our study was 53 protein-ligand complexes obtained from the RCSB Protein data bank (PDB) (Burley et al., 2018). The specific PDB access codes for these structures can be found in Table S1. Initial preparation was conducted on the dataset to ensure consistency and to address potential variations resulting from different protocol implementations in different software. While this section covers general preparation steps, software-specific tasks, such as format conversions, are elaborated upon in their respective subsections. A unique case within our dataset was system 6CNK (Walsh et al., 2018), a heteromeric (α4)3(β2)2 pentamer featuring binding sites occupied at two distinct interfaces (α4/α4 and α4/β2). To accurately represent the complexity of
this system, it was divided into two separate entries, denoted as 6CNKAB and 6CNKBC, thereby increasing the total number of systems studied to 54. During our statistical analysis, we considered different categories reflecting the nature of the protein targets and ligands. We defined two primary groups: the “neonicotinoids-like” group, which includes compounds directly derived from neonicotinoids, sulfoximines, and butenolides; and the “drug-like” group, encompassing small molecules that do not relate to the neonicotinoid category such as lobeline and cytisine. Furthermore, some AChBPs were reclassified into the nAChR group because of humanization or insectification of their binding sites. The detailed classifications and associated PDB codes are given in Table S1 and Fig. S1.

### 2.1.1. Protein and receptors preparation

The protein crystallographic structures were prepared using HTMD system preparation workflow (Doerr et al., 2016). The protonation states were computed via ProPka 3 (Olsson et al., 2011), after which bond orders, formal charges, and hydrogens were systematically added to each system. Special attention was paid to the tautomeric states of histidine residues, as well as to the rotamers of asparagine and glutamine residues as adjustments were made to ensure an optimized hydrogen bonding network. While cocrystallized water molecules were generally removed, exceptions were made for those within a 3 Å radius of the bound ligands. Specific cases, notably 5AFH (Spurny et al., 2015) and 7KOQ (Noviello et al., 2021), presented missing disulfide bonds in the loop C region due to crystallization artifacts, and these were rectified manually. Each structure underwent energy minimization using the CHARMM36m force field (Huang et al., 2016) and NAMD 2.14 (Kalé et al., 1999; Phillips et al., 2005, 2020). During this process, sidechain atoms were permitted unrestricted movement, whereas backbone atoms were constrained with a harmonic restraint of 5 kcal/mol/Å².

### 2.1.2. Ligands preparation

The compounds addressed in this study were initially identified and extracted from the aforementioned 53 protein-ligand complexes. They were then accurately reconstructed from SMILES string formulas utilizing RDKit's EmbedMolecule function, adhering to default parameters (Landrum, 2021). Subsequent minimization of these structures was achieved using the MMFF94 (Halgren, 1996) implementation in RDKit. For each compound, tautomers and protonation states were generated at pH 7.4 using Epik v5.5 (Schrödinger release, 2021). Once the correct protonation states and tautomer conformations identified, the ligands geometries were optimized at the second-order Møller-Plesset perturbation theory (MP2) level of theory in combination with the 6-31G(d) basis set using Gaussian 16 (Frisch et al., 2016).

## 2.2. Molecular docking

For our molecular docking analyses, we employed several established software tools: Ledock v1.0, ADFR suite 1.1, Autodock Vina 1.2 (referred to as Vina henceforth), and the CDOCKER pycharmm implementation version 1.0 (referred to as CDOCKER henceforth). Both semi-flexible and flexible docking protocols were investigated across all software, except for Ledock, which does not support flexible docking. In the context of this study, semi-flexible docking refers to a protocol where the ligand is allowed to move freely, but the protein receptor remains rigid. In contrast, flexible docking permits full adaptability of both entities. For semi-flexible docking approaches, the binding site region was delineated within a cube centered on the endogenous ligand, with each side measuring 6 Å. For the flexible docking protocol, this cube was enlarged to 8 Å to accommodate greater potential movement. However, only the amino acids from the initial 6 Å cube were chosen for flexibility. Post-docking, the results were visualized using Pymol 2.5.0.

### 2.2.1. Ledock

The prepared structures were adjusted to correspond to CHARMM36m protein force field naming scheme (Huang et al., 2016). A 6 Å cube, centered on the ligand, was established as the docking box to encompass the binding site, with coordinates generated using the GetBox-PyMOL-Plugin (Xiao, 2018). We set the output to yield 10k structures and used 0.5 Å as internal clustering threshold. Ligand structures were converted to the necessary mol2 format using OpenBabel version 3.1.0 (O'Boyle et al., 2011; The open babel package, 2020).

### 2.2.2. ADFR

Both ligand and receptor structures, initially processed by the HTMD module, were subsequently refined using ADFR suite v1.1rc1 to ensure they were aptly prepared for ADFR. For the receptor, the *prepare_receptor* script was employed, generating the desired pdbqt file format. During this process, Gasteiger partial charges were computed, and non-polar hydrogens, along with dummy atoms symbolizing lone pairs, were discarded. The partial charges associated with these elements were then merged with their associated heavy atom. In parallel, ligands were processed using the *prepare_ligand* script, with all rotatable bonds granted the freedom to rotate.

For semi-flexible docking, we used Autogrid v4.2.7 to generate affinity maps for every atom type present in the ligand set, desolvation maps, and electrostatic maps. These maps were based on docking center coordinates with specific grid settings: 0.375 Å spacing and 0.5 smoothing value. We utilized the Lamarckian Genetic Algorithm (LGA) for conformation generation, setting the genetic algorithm (GA) to run 10k times with a 2.5M energy evaluation cap per iteration. For flexible docking, parameters mirrored those of semi-flexible docking, with one significant deviation: the specification of flexible amino acids within the binding site. Notably, cysteines constituting the disulfide bridge in the loop C region were omitted due to the absence of rotamer data for them (Dunbrack and Karplus, 1994).

### 2.2.3. Autodock Vina

For our study, ligands underwent a preparation process using ADFR suite v1.1rc1, mirroring the methods applied in the **ADFR** section. In the context of Vina, we explored three force fields: ad4 (AutoDock 4's implementation within Vina), Vina, and Vinardo. For the ad4 force field, receptor preparation paralleled the **ADFR** approach, with grid maps generated via Autogrid 4.2.7. Conversely, when employing Vina and Vinardo force fields, no pre-computed maps were necessary as these are generated in real-time during the docking process. Common parameters across all three force fields included an exhaustiveness of 512, a maximum output of 20 binding modes, and a maximum energy difference of 3 kcal/mol between modes. The default parameters were maintained for the respective scoring function weights and terms.

For the flexible docking strategy, receptors were processed into two sets of files: one featuring the rigid components and another containing only the amino acids considered for flexibility. For ad4 force field, a method akin to the **ADFR** section for grid map generation was followed. In contrast, Vina and Vinardo force fields computed these maps on-the-fly. It's important to note that cysteines forming the disulfide bridge in the loop C region were omitted due to the absence of relevant rotamer data (Dunbrack and Karplus, 1994).

### 2.2.4. CDOCKER

Using CDOCKER, which is based on the CHARMM molecular dynamics framework, our ligands were processed with CGenFF v2.5 first (Vanommeslaeghe and MacKerell, 2012; Vanommeslaeghe et al., 2012). This processing generated CHARMM compatible parameters, which were saved as stream files for each ligand. Both ligands and receptors were subjected to a naming scheme matching to ensure compatibility
with CGenFF force field v4.6 (Vanommeslaeghe et al., 2010) and CHARMM36m protein force field (Huang et al., 2016). CHARMM version c47b2 (Brooks et al., 2009) facilitated the generation of the essential input files: .pdb for coordinates, .psf for topology, and .rtf and .prm for parameters.

For the semi-flexible docking phase, OpenBabel 3.1.0 was employed to produce 500 starting conformations. These conformations were then positioned in 25 distinct locations within the predefined binding site, yielding a total of 125k poses for each ligand. CDOCKER minimization was utilized to refine these poses, ensuring their optimal fit within the binding site. Specific parameters were established, including a dielectric constant set to 3.0, “high” exhaustiveness, and the use of the hydrogen bond grid potential. Default settings were retained for parameters associated with the soft, hard, and native grids, inclusive of their cutoffs. The flexible docking approach mirrored this process but utilized a reduced number of starting positions: 20 initial conformers, 25 placements, and 50 iterations for the genetic algorithm.

## 2.3. Statistical analysis

To evaluate the linear association between the optimal root-mean-square deviation (RMSD) from reference crystallographic structures and the docking scores, Pearson correlation coefficients for each dataset were first computed. Due to the heterogeneity in dataset sizes, a weighted Fisher-z transformation (Fisher, 1992) was employed to normalize the coefficients, thereby allowing for the aggregation of data from all datasets. This normalization was necessary as the reliability of Pearson correlation coefficients decreases with the size of the dataset, particularly for those with fewer than three data points. The Fisher-z transformation of a Pearson coefficient *r* is given by:

$$z = \frac{1}{2} \ln \left( \frac{1+r}{1-r} \right) \qquad (1)$$

where *r* is the Pearson correlation coefficient and *z* is the Fisher-z value. Each *z* value was weighted by the sample size *N* of its dataset minus three (*N* − 3), considering the standard deviation of the Fisher-z distribution is $\frac{1}{\sqrt{N-3}}$:

$$z_w = z \times (N - 3) \qquad (2)$$

The aggregation of these weighted values yields a mean weighted Fisher-z value $\bar{Z}_w$, which, upon applying the inverse Fisher-z transformation, returns a mean Pearson correlation coefficient $\bar{r}$ reflective of the combined datasets:

$$\bar{r} = \frac{e^{2\bar{Z}_w} - 1}{e^{2\bar{Z}_w} + 1} \qquad (3)$$

The standard error for $\bar{Z}_w$ was determined by the reciprocal of the square root of the sum of the weights:

$$SE_{\bar{Z}_w} = \left( \sum_{i=1}^{k} (N_i - 3) \right)^{-\frac{1}{2}} \qquad (4)$$

where *k* is the number of datasets, and $N_i$ is the sample size for the *i*<sub>th</sub> dataset.

Thus, the z-score ($Z_{score}$) for $\bar{Z}_w$ is given by:

$$Z_{score} = \frac{\bar{Z}_w}{SE_{\bar{Z}_w}} \qquad (5)$$

Lastly, the p-values were subsequently derived from the cumulative distribution function ($\Phi$) of the standard normal distribution, reflecting the probability of observing a z-score as extreme as or more extreme than what was calculated under the null hypothesis of no correlation. This was performed as a two-tailed test, given by:

$$p = 2 \times (1 - \Phi(|Z_{score}|)) \qquad (6)$$

To investigate the data's non-linear and non-parametric correlation characteristics, we computed Spearman's rank correlation coefficient and Kendall's $\tau$ coefficient. These rank-based statistics are resilient to the effects of outliers and do not presuppose normality of data distribution, rendering the Fisher-z transformation inapplicable. Instead, we implemented a bootstrap resampling method to ascertain the average correlation coefficients and their corresponding confidence intervals. By performing 10k bootstrap iterations and targeting a 99 % confidence level, we obtained a distribution for each coefficient that does not rely on theoretical distribution models.

# 3. Results and discussion

## 3.1. Semi-flexible docking: assessment against crystal structures

In our study, the investigation into the reproducibility of docking poses, as generated by a diverse array of docking protocols across multiple tools, was a crucial component. This analysis was conducted against crystallographic reference poses, offering a window into each tool's capability to accurately predict the binding modes of compounds. This aspect was especially pertinent for insecticides, known for their unique binding characteristics. The RMSD served as our primary metric for this evaluation, with docked poses inspected against thresholds of 0.5 Å, 1.0 Å, and 1.5 Å. These thresholds were deliberately chosen to reflect different levels of rigor in assessing pose accuracy. Our adoption of the more stringent lower threshold of 0.5 Å, which deviates from the conventional benchmarks commonly cited in existing literature (Ravindranath et al., 2015; Eberhardt et al., 2021; Wang et al., 2016), stemmed from an explicit requirement for accurate exactness. This threshold was essential for accurately describing the specific binding modes of insecticides, a factor of high importance given their complex interaction dynamics and potential environmental and biological implications and is used as main accepted threshold for the upcoming discussion.

### 3.1.1. Accuracy assessment

A central aspect of our analysis involved a detailed examination of each descriptor characterizing optimal docking results within the semi-flexible docking approach. For each protocol and software, we assessed the ability to yield accurate results against crystallographic references, independent of pose ranking and corresponding scores. Consequently, from the array of results returned for each system, we singled out only the best pose, calculating its RMSD in comparison to the reference structure. The results of this process are illustrated in Fig. 1, and an accuracy percentage derived from these results is presented in Fig. 2. We classified the accuracy measurements into nine distinct categories, reflecting the nature of the target (AChBP or nAChR), the type of ligand (drug-like or neonicotinoid-like), various combinations of these two factors, and an overarching total category encompassing all data.

An initial review of the data reveals significant disparities dependent on the chosen RMSD threshold. At the stringent 0.5 Å threshold, all software tools demonstrated relatively low accuracy. Within the total category, a descending order of accuracy was observed as follows: Ledock, CDOCKER, Vina_vinardo, Vina_vina, Vina_ad4, and ADFR, with respective accuracies of 31 %, 26 %, 17 %, 17 %, 11 %, and 9 %. Analyzing the other categories, we noted that Ledock and CDOCKER exhibited a similar trend, demonstrating relative effectiveness, especially with AChBP targets and neonicotinoid ligands. Conversely, ADFR, Vina_vina, Vina_ad4, and Vina_vinardo displayed less distinction between categories, with the exception of Vina_vinardo, which showed a preference for nAChR targets, though no significant differentiation was observed between ligand types. Notably, in the nAChR neonicotinoid category, almost all software tools faced significant challenges at this threshold, although the statistical significance of this observation is moderated by the category's size.

At the more lenient 1.5 Å threshold, an overall excellent accuracy
<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Scatter plot illustrating the best RMSD values of docked poses against the respective system names, segregated by software for the semi-flexible docking method. The colored sections categorizes the systems based on the nature of the ligand (drug-like or neonicotinoid-like) and the type of receptor (AChBP or nAChR). Horizontal dashed lines represent RMSD thresholds that defines different levels of docking accuracy.</figcaption></figure>

was observed across all software tools. The accuracy ranking for the total category in descending order is as follows: Ledock, Vina_ad4, Vina_vinardo, Vina_vina, CDOCKER, and ADFR, with accuracies of 89 %, 81 %, 65 %, 65 %, 57 %, and 52 %, respectively. In broader ligand-based categories such as drug-like and neonicotinoid-like, the performance variances of the software tools were similarly revealing. With the exception of ADFR at this threshold, all tools showed marked proficiency, suggesting that the nature of the ligand exerts a significant influence on docking accuracy. In contrast, for protein targets, a higher accuracy was observed with nAChR across all software. In summarizing the overall performance across these categories, Ledock consistently exhibited the highest accuracy, irrespective of the chosen threshold. This consistent performance underscores Ledock's utility in diverse docking scenarios, affirming its robustness and versatility in accurately predicting docking poses.

### 3.1.2. Precision assessment: correlation between RMSD values and docking scores

In the continuation of our comprehensive analysis of docking software performances on our datasets, we turned our attention to a crucial aspect complementary to accuracy: precision. The standout accuracy of Ledock among the tested software prompted an in-depth evaluation of the effectiveness of scoring functions in accurately ranking docked poses. This data is given in Fig. 3, where we analyzed the correlation between the RMSD values of docked poses and their corresponding scores. Our methodology involved the use of the weighted average Pearson correlation coefficient, computed using Fisher's Z transformation for linear correlation characterization, alongside the bootstrapped average Spearman and Kendall τ correlation coefficients for monotonic relationship assessment. These statistical aggregation techniques were employed to synthesize the data across each software and protocol, while mitigating the influence of their varying sample size outputs. Analyses were conducted for both the complete output poses and the top twenty poses subsets, ensuring a comprehensive measure of each tool's precision capabilities.

The findings revealed that the weighted average Pearson correlation coefficients for all software tools were significantly lower than the generally accepted values indicative of a strong correlation. Interestingly, Ledock, despite its high accuracy, emerged as the least precise. Its effective correlation between poses RMSD and scores was particularly low, at 0.055 for the top twenty poses subset and 0.117 for the complete

output poses subset, indicating a markedly ineffective linear correlation. The corresponding p-values of $1.25 \times 10^{-21}$ and $1.77 \times 10^{-42}$ for both respective subsets reinforce the absence of a meaningful linear correlation. Similarly, while ADFR and Vina_ad4 showed slightly higher weighted average Pearson coefficients than Ledock, they did not achieve satisfactory levels of precision. Vina_vina and vina_vinardo, despite having the highest weighted average Pearson coefficients, still fell short, highlighting a general limitation in their precision. CDOCKER, aligning with the other tools, exhibited Pearson coefficients well below the threshold for strong correlation, especially in the top twenty poses subset, thus indicating limited precision.

Further insights were gleaned from the non-parametric rank correlations, measured by the bootstrapped Spearman and Kendall τ coefficients. These coefficients revealed a weak rank-order relationship for Ledock and CDOCKER, with broad confidence intervals indicating significant uncertainty in their correlation strength, irrespective of the selected subset. ADFR and Vina_ad4 fared slightly better in these measures, yet still pointed to a weak rank-order relationship, albeit with more consistent correlation patterns. Vina_vina and Vina_vinardo, despite registering the highest coefficients, exhibited a similar trend, denoting a weak rank-order relationship and a lack of consistent reliability in their scoring functions' ability to rank docked poses. The implications of these weak non-parametric correlations are substantial. They suggest that the scoring functions of these tools are not consistently ranking docked poses in a way that correlates with their actual RMSD values. This inconsistency is of particular concern in docking studies where precision is crucial. Top-ranked poses by these scoring functions may not necessarily be the closest to the reference structure in terms of RMSD, potentially leading to misleading conclusions. Therefore, when juxtaposing these precision findings with the accuracy results from the previous section, an intriguing paradox emerges. The most accurate tool, Ledock, ironically turns out to be the most random in terms of score-to-RMSD correlation, whereas Vina_vinardo, ranking third in accuracy, appears to be the most consistent. These observations underscore the complexity and challenges in selecting the appropriate software and scoring function for docking studies involving neonicotinoid-like compounds.

### 3.1.3. Docking score against experimental binding affinities

In our comprehensive exploration of docking software capabilities, we honed in on a crucial aspect of this study: elucidating the relationship
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Bar plot depicting the percentage of valid docked poses below specific RMSD thresholds, categorized by docking software and the class of receptor or compound for semi-flexible docking method.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Lollipop plot displaying the aggregated correlation coefficients between RMSD values and docking scores per software for semi-flexible docking method. For each software, the plot presents the weighted average Pearson correlation coefficient after Fisher-z transformation, as well as the bootstrap resampling averages for Spearman and Kendall τ coefficients. Colored markers denote the correlation for all poses, while the ‘X’ markers indicate the correlation for the top twenty poses by score. In the Pearson plot, p-values are displayed adjacent to the markers, presented in the format (top twenty subset, complete output poses subset). Diamond markers represent the confidence intervals for the complete output poses subset, and black bars denote the confidence intervals for the top twenty poses subset.</figcaption></figure>

between docking scores and experimental *K*<sub>d</sub> values. Due to the absence of experimental data for all 54 systems, we focused on a select subset of 15 systems for which data were available (from the PDBbind database (Wang et al., 2005)). This selection is detailed in the Table S1. Fig. S4 illustrates the docking scores versus experimental *K*<sub>d</sub>, categorized into three groups: All data, neonicotinoid-like, and drug-like. Employing

correlation coefficients (Pearson, Spearman, and Kendall τ) for each software tool across these categories, the results are presented in Fig. 4, depicted as a lollipop plot with statistical analysis.

Ledock displayed a moderately low but consistent predictive capability in the all data subset, with Pearson, Spearman, and Kendall τ correlations around 0.4, each statistically significant (p-values: 0.13).
<figure><img src="image_8.png" /><figcaption>Fig. 4. Lollipop plots depicting the correlation between docking scores of the best poses and experimental ($K_d$) values for various docking software using the semi-flexible method. The left panel shows the Pearson correlation coefficients, the middle panel shows the Spearman correlation coefficients, and the right panel shows the Kendall $\tau$ correlation coefficients. Colored markers represent the entire dataset (both drug-like and neonicotinoid-like classes), while the grey cross markers denote the neonicotinoid-like class and the beige plus markers denote the drug-like class specifically. The annotations provide the p-value for all dataset (A), for drug-like subset (D) and the neonicotinoid-like subset (N) for each software.</figcaption></figure>

However, in the drug-like subset, a strikingly high Pearson correlation of 0.940 (p-value: 0.001) was observed, a stark contrast to the lower Spearman (0.690, p-value: 0.058) and Kendall $\tau$ (0.571, p-value: 0.061) correlations. While the significant p-values indicate a robust linear relationship, the discrepancy with rank-based correlations suggests potential biases, possibly influenced by a small number of high-performing compounds. Conversely, the neonicotinoid-like subset exhibited

negative correlations across all metrics, with the p-values suggesting these relationships are not statistically significant, reflecting the unpredictable nature of these interactions. ADFR showed a pattern similar to Ledock in the all data subset, with a moderate Pearson correlation of 0.539 (p-value: 0.038) that was not echoed in the lower Spearman (0.296, p-value: 0.283) and Kendall $\tau$ (0.181, p-value: 0.379) correlations. The significant Pearson p-value, in contrast with the insignificant

<figure><img src="image_9.png" /><figcaption>Fig. 5. Bar plot depicting the percentage of valid docked poses below specific RMSD thresholds, categorized by docking software and the class of receptor or compound for flexible docking method.</figcaption></figure>
rank-based correlations, points to a non-linear relationship possibly influenced by non-uniform distributions. The Drug-like subset, again, showed a robust Pearson correlation of 0.958 (p-value: 0.001).

Vina_ad4, Vina_vina, and Vina_vinardo exhibited moderate to low Pearson correlations in the all data subset, with variably lower and less significant Spearman and Kendall τ correlations. This pattern indicated nuanced relationships between docking scores and binding affinities, not fully captured by linear correlations. In the Drug-like subset, despite the significant p-values indicating strong linear correlations, caution is warranted due to the lower rank-based correlations, suggesting a potential linear bias.

CDOCKER consistently showed the lowest performance, with negative Pearson (− 0.348, p-value: 0.203), Spearman (− 0.452, p-value: 0.196), and Kendall τ (− 0.219, p-value: 0.275) correlations in the all data subset. Contrary to other software, the drug-like subset was poorly correlated, as indicated by negative coefficients and statistically insignificant p-values. The neonicotinoid-like subset consistently revealed lower or negative correlations across all metrics for all software and protocols, with the p-values often denoting a lack of statistical significance. This underscores the complexity and challenges in establishing correlations between docking results and experimental data for this specific class of compounds.

## 3.2. The evaluation of flexible docking

The docking protocols discussed earlier primarily utilized the semi-flexible docking approach, a mainstream method in molecular simulations for such systems. Expanding on this, we explored the flexible docking approach, where both ligand and receptor exhibit flexibility. Theoretically, flexible docking, due to its potential to better capture the induced-fit effect, is expected to yield superior results compared to semi-flexible docking (Kitchen et al., 2004).

Data from Figs. S2, S3 and 5 provide a comprehensive overview of the accuracy achieved by all flexible docking software and protocols, excluding Ledock which lacks a flexible docking option. At the stringent threshold of 0.5 Å, a notable lack of accuracy is evident in the results. Vina_Vina emerges as the most accurate in the total category with 7 % accuracy, followed by ADFR, Vina_ad4, and Vina_vinardo with 4 %, 6 %, and 2 % respectively. Remarkably, CDOCKER fails to achieve successful docking for all systems, leading to a 0 % accuracy rate. When the threshold is relaxed to 1.5 Å, there is a marked improvement in accuracy. ADFR and Vina_vina achieve similar accuracies of 46 % and 44 %, while Vina_ad4 and Vina_vinardo reach 33 % and 19 % respectively. Consistent with earlier observations, CDOCKER maintains a 0 % accuracy rate. The trends across categories remain inconclusive, although a pattern similar to that in the semi-flexible docking accuracy assessment is observed, where the neonicotinoid-like category appears to yield more successful results than the drug-like one for the more accurate

approaches. Precision wise, as presented in Fig. 6, shows a decrease with the flexible docking approach. The weighted average Pearson correlation coefficients and their p-values, along with bootstrapped averaged Spearman and Kendall τ values, reflect this. ADFR displays modest Pearson coefficients (0.244 for complete output poses subset, 0.119 for top twenty subset), suggesting a poor linear correlation. Vina_ad4, with very low Pearson coefficients (0.023 for complete output poses subset, 0.005 for top twenty subset) with low p-values (0.0002), indicates a weak, albeit statistically significant, lack of linear correlation. Vina_vina, Vina_vinardo, and CDOCKER exhibit very low Pearson coefficients, reinforcing the absence of linear correlation between RMSD and scores. The rank-based coefficients mirror this trend, with values around zero for Spearman and Kendall τ for all software, except for Vina_ad4 and CDOCKER, where the values are negative, hinting at a complete absence of correlation. The wide intervals for each software reinforce this observation for both subsets. The following observation is also reflected in the correlation between docking scores and experimental $K_d$ values. The corresponding data are given in Figs. S5 and 7. ADFR exhibited a low Pearson correlation on the total dataset, with a coefficient around 0.499, which was statistically significant (p-value: 0.058). However, the drug-like subset showed an exceptionally high Pearson correlation of 0.970 (p-value: 0.001), contrasting with lower Spearman (0.405, p-value: 0.320) and Kendall τ (0.214, p-value: 0.548) correlations. This stark difference implies a strong linear relationship for drug-like compounds, although the lower rank-based correlations suggest potential biases or non-linearity in the data. In contrast, the neonicotinoid-like subset displayed lower correlations across all metrics, reflecting the challenging and unpredictable nature of these interactions, as evidenced by the Pearson correlation of 0.271 (p-value: 0.556). Other software tools, including Vina_ad4, Vina_vina, Vina_vinardo, demonstrated varying degrees of correlation in the total data subset, with a general trend of low Pearson correlations and even lower Spearman and Kendall τ correlations. CDOCKER notably showed consistently lower performance across all datasets, with negative Pearson, Spearman, and Kendall τ correlations in the total dataset. This pattern was also observed in the drug-like subset, indicating poor correlation and statistical insignificance in p-values. The neonicotinoid-like subset further highlighted the challenges in establishing meaningful correlations for this class of compounds across all software and protocols, often exhibiting negative correlations and statistically insignificant p-values, underscoring the inherent complexity in modeling these correlations.

Despite the initial concerns raised by these seemingly counterintuitive findings, it is important to clarify that these results are not attributable to insufficient sampling of the conformational space. Rather, they reflect the inherent limitations in the accuracy and precision of flexible docking. Throughout this study, we meticulously selected a set of parameters to ensure an extensive exploration of the conformational space for each ligand-protein complex. Specifically, considering the largest

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Lollipop plot displaying the aggregated correlation coefficients between RMSD values and docking scores per software for flexible docking method. For each software, the plot presents the weighted average Pearson correlation coefficient after Fisher-z transformation, as well as the bootstrap resampling averages for Spearman and Kendall τ coefficients. Colored markers denote the correlation for complete output poses subset, while the 'X' markers indicate the correlation for the top twenty poses by score. In the Pearson plot, p-values are displayed adjacent to the markers, presented in the format (top twenty dataset, complete output poses subset). Diamond markers represent the confidence intervals for the complete output poses subset, and black bars denote the confidence intervals for the top twenty poses dataset.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Lollipop plots depicting the correlation between docking scores of the best poses and experimental ($K_d$) values for various docking software using the flexible method. The left panel shows the Pearson correlation coefficients, the middle panel shows the Spearman correlation coefficients, and the right panel shows the Kendall $\tau$ correlation coefficients. Colored markers represent the entire dataset (both drug-like and neonicotinoid-like classes), while the grey cross markers denote the neonicotinoid-like class and the beige plus markers denote the drug-like class specifically. The annotations provide the p-value for all dataset (A), for drug-like subset (D) and the neonicotinoid-like subset (N) for each software.</figcaption></figure>

ligand in our dataset, which features 13 rotatable bonds, and combining this with an average of 8 flexible residues per binding site–each possessing approximately 2 rotatable bonds–the total number of rotatable bonds to explore reaches 29. Given this complexity, the default parameters, typically deemed adequate for systems with up to 10 rotatable bonds for Vina (Eberhardt et al., 2021; C. of Computational Structural Biology (CCSB) Scripps Research, 2021) and up to 20 for ADFR (Ravindranath et al., 2015; C. of Computational Structural Biology (CCSB), 2019), were significantly adjusted. For Vina-based protocols, exhaustiveness was increased from 8 to 512, coupled with a maximum evaluation limit of 2.5M. For ADFR, the initial genetic algorithm generation parameters were enhanced from 50 to 10k, and the maximum evaluation per generation was increased from 20k to 2.5M, enabling a total of 25 billion dockings. For CDOCKER's pyCHARMM implementation, the optimal number of starting conformations that can fit within the binding site was identified and set to 25k, as generating more than this results in initial conformations that extend beyond the confines of the binding region. In all instances, the maximum allowed number of evaluations was never reached, ensuring comprehensive coverage of the potential conformational space. This extensive exploration reinforces the conclusion that observations related to flexible docking performance are robust and not merely artifacts of inadequate sampling.

In summary, while the flexible docking approach increases complexity, it fails to enhance accuracy and precision. This suggests that the added computational intricacy may not necessarily translate into improved predictive performance for our compounds, a finding that resonates with similar studies on other environmental chemicals (Xue et al., 2022).

### 3.3. Computational cost

The computational cost analysis, comparing flexible and semi-flexible docking approaches, offers vital insights into the performance and efficiency of various docking software in environmental chemical studies. For this study, we used a node with a 2.40 GHz Intel Xeon Silver 4210R CPU with 20 cores for calculation involving Ledock, ADFR, Vina_ad4, Vina_vina, Vina_vinardo software for both semi-flexible and flexible docking. Since CDOCKER is running on GPU, a different node with a 2.50 GHz Intel Xeon Gold 6248 CPU with 20 cores and a Tesla V100 GPU was used for both docking methods involving this software.

The data, presented in Fig. S6, provides an understanding of the run times associated with each software under both docking conditions. In the semi-flexible docking subset, most software tools, notably Vina_vina and Vina_vinardo, demonstrated efficient performance with relatively lower run times (under one hour). Ledock and ADFR are notably slower as the average run time is over 10 h per ligand. However, we can highlight that the execution time is highly system dependent.

In contrast, the flexible docking approach, as detailed in the second subset, showed a substantial increase in run times across all software tools. Table S3 further clarifies this by showing the ratio of run times

between flexible and semi-flexible approaches for each software, revealing stark increases in computational time. Vina_ad4, Vina_vina and Vina_vinardo for example, showed dramatic increases, with ratios of 108.86, 92.27 and 77.43 respectively. ADFR also exhibited notable increases in computational requirements, with a ratio of 14.01. CDOCKER seems to be the more consistent between both protocols as the increase in run time is only of 1.82 folds. These numbers accentuate the significant rise in computational load associated with flexible docking, which becomes particularly relevant and ironic when considering the absence of improvements, in predictive performance.

### 4. Conclusions

In this study, we embarked on a detailed exploration of various prominent docking software, employing both semi-flexible and flexible docking paradigms. Our objective was to critically assess their ability to reproduce experimental results, particularly crystallographic data, using robust statistical approaches. This analysis has provided valuable insights into the effectiveness of these tools in simulating molecular interactions, especially in the context of insecticides. In the realm of semi-flexible docking, Ledock emerged as the most accurate software, while Autodock Vina, employing the Vinardo scoring function, proved to be the most reliable. However, when considering both accuracy and reliability, the picture remained unclear, indicating that no single software could be deemed superior in all aspects. This underscores the complexities involved in computational docking and the need for careful consideration of the specific requirements of each study. A pivotal aspect of our analysis focused on the correlation between docking scores and experimental $K_d$ values. Interestingly, none of the software tools demonstrated a clear correlation for neonicotinoid-like compounds. In contrast, drug-like compounds exhibited excellent correlation, as per the available data. This stark difference highlights a significant bias in the scoring functions of existing software, which appear to be tailored predominantly towards drug-like compounds. This finding is crucial, as it implies that these tools may not be suitable for environmental chemicals, specifically insecticides, that do not conform to the criteria of therapeutic drugs. Furthermore, our comparison between semi-flexible and flexible docking paradigms revealed that flexible docking, despite its increased computational complexity, failed to yield more accurate results. This was compounded by a significant escalation in computational costs, rendering this approach virtually ineffective for the class of compounds under study. The increased complexity and resource demand associated with flexible docking did not translate into better predictive accuracy or reliability, particularly for neonicotinoid-like compounds. Therefore, based on our findings, it is imperative for the scientific community to continue the characterization and evaluation of available docking tools. More importantly, there is a clear need to consider the development of specialized tools tailored for chemicals that do not align with the characteristics of conventional therapeutic drugs. Such advancements are crucial for enhancing our understanding and
predictive capabilities in the realm of environmental chemoinformatics, ultimately contributing to more informed research in this critical area.

## Declaration of Competing Interest

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

## Funding

This research was funded by the French Agency for Food, Environmental and Occupational Health and Safety (ANSES), Grant no. 2020/01/214project EVALSUB.

## CRediT authorship contribution statement

**Zakaria Bouchouireb:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.
**Damien Olivier-Jimenez:** Investigation. **Titouan Jaunet-Lahary:** Investigation. **Steve H. Thany:** Writing – review & editing, Funding acquisition. **Jean-Yves Le Questel:** Writing – review & editing, Funding acquisition.

## Data availability

Data will be made available on request.

## Acknowledgments

This work was performed using HPC resources from (i) CCIPL (Centre de Calcul Intensif des Pays de Loire), (ii) GenOuest bioinformatics core facility (https://www.genouest.org) and (iii) GENCI-IDIRS (Grant 2023-0140714125). They are acknowledged for providing the computing infrastructure and computing time.

## Appendix A. Supporting information

Supplementary data associated with this article can be found in the online version at doi:10.1016/j.ecoenv.2024.116582.

## References

Alamididine, Z., Selvam, B., Graton, J., Laurent, A.D., Landagaray, E., Lebreton, J., Mathé-Allainmat, M., Thany, S.H., Questel, J.-Y.L., 2019. Binding of sulfoxaflor to *Aplysia californica* achnb: computational insights from multiscale approaches. J. Chem. Inf. Model. 59 (9), 3755-3769. https://doi.org/10.1021/acs.jcim.9b00272.

Bael, J.B., Holloway, G.A., 2010. New substructure filters for removal of pan assay interference compounds (PAINS) from screening libraries and for their exclusion in bioassays. J. Med. Chem. 53 (7), 2719-2740. https://doi.org/10.1021/jm901137j.

Brooks, B.R., Brooks, C.L., Mackerell, A.D., Nilsson, L., Petrella, R.J., Roux, B., Won, Y., Archontis, G., Bartels, C., Boresch, S., Caflisch, A., Caves, L., Cui, Q., Dinner, A.R., Feig, M., Fischer, S., Gao, J., Hodoscek, M., Im, W., Kuczera, K., Lazaridis, T., Ma, J., Ovchinnikov, V., Paci, E., Pastor, R.W., Post, C.B., Pu, J.Z., Schaefer, M., Tidor, B., Venable, R.M., Woodcock, H.L., Wu, X., Yang, W., York, D.M., Karplus, M., 2009. CHARMM: the biomolecular simulation program. J. Comput. Chem. 30 (10), 1545-1614. https://doi.org/10.1002/jcc.21287.

Buckner, J., Liu, X., Chakravorty, A., Wu, Y., Cervantes, L.F., Lai, T.T., Brooks, C.L., 2023. pyCHARMM: embedding CHARMM functionality in a python framework. J. Chem. Theory Comput. 19 (12), 3752-3762. https://doi.org/10.1021/acs.jctc.3c00364.

Burley, S.K., Berman, H.M., Bhikadiya, C., Bi, C., Chen, L., Costanzo, L.D., Christie, C., Dalenberg, K., Duarte, J.M., Dutta, S., Feng, Z., Ghosh, S., Goodsell, D.S., Green, R. K., Guranović, V., Guzenko, D., Hudson, B.P., Kalro, T., Liang, Y., Lowe, R., Namkoong, H., Peisach, E., Periskova I., Prlić, A., Randle, C., Rose, A., Rose, P., Sala, R., Sekharan, M., Shao, C., Tan, L., Tao, Y.-P., Valasatava, Y., Voigt, M., Westbrook, J., Woo, J., Yang, H., Young, J., Zhuravleva, M., Zardecki, C., 2018. RCSB protein data bank: biological macromolecular structures enabling research and education in fundamental biology, biomedicine, biotechnology and energy. Nucleic Acids Res. 47 (D1), D464-D474. https://doi.org/10.1093/nar/gky1004.

C. of Computational Structural Biology (CCSB) Scripps Research, 2021. Revision d13f9e93, Autodock vina: Molecular docking program: Documentation. https://autodock-vina.readthedocs.io/en/latest/index.html), (Accessed 06 March 2024).

C. of Computational Structural Biology (CCSB), 2019. Autodockfdr documentation, https://ccsb.scripps.edu/adfr/documentation/, (Accessed 06 March 2024).

Cartereau, A., Taillebois, E., Questel, J.-Y.L., Thany, S.H., 2021. Mode of action of neonicotinoid insecticides imidacloprid and thiacloprid to the cochwalla pame *α*7 nicotinic acetylcholine receptor. Int. J. Mol. Sci. 22 (18), 9880. https://doi.org/10.3390/ijms22189880.

Casida, J.E., Durkin, K.A., 2013. Neuroactive insecticides: targets, selectivity, resistance, and secondary effects. Annu. Rev. Entomol. 58 (1), 99-117. https://doi.org/10.1146/annurev-ento-120811-153645.

Chaput, L., Martinez-Sanz, J., Saettel, M., Mouawad, L., 2016. Benchmark of four popular virtual screening programs: construction of the active/decoy dataset remains a major determinant of measured performance. J. Cheminf. 8 (1) https://doi.org/10.1186/s13321-016-0167-x.

Chen, L., Cruz, A., Ramsey, S., Dickson, C.J., Duca, J.S., Hornak, V., Koes, D.R., Kurtzman, T., 2019. Hidden bias in the dude dataset leads to misleading performance of deep learning in structure-based virtual screening. PLoS One 14 (8), e0220113. https://doi.org/10.1371/journal.pone.0220113.

Chen, D., Zhang, Y., Lv, B., Liu, Z., Han, J., Li, J., Zhao, Y., Wu, Y., 2020. Dietary exposure to neonicotinoid insecticides and health risks in the chinese general population through two consecutive total diet studies. Environ. Int. 135, 105399 https://doi.org/10.1016/j.envint.2019.105399.

Cheng, T., Li, Q., Zhou, Z., Wang, Y., Bryant, S.H., 2012. Structure-based virtual screening for drug discovery: a problem-centric review. AAPS J. 14 (1), 133-141. https://doi.org/10.1208/s12248-012-9322-0.

Coleman, R.G., Carchia, M., Sterling, T., Irwin, J.J., Shoichet, B.K., 2013. Ligand pose and orientational sampling in molecular docking. PLoS One 8 (10), e75992. https://doi.org/10.1371/journal.pone.0075992.

Cross, J.B., Thompson, D.C., Rai, B.K., Baber, J.C., Fan, K.Y., Hu, Y., Humblet, C., 2009. Comparison of several molecular docking programs: pose prediction and virtual screening accuracy. J. Chem. Inf. Model. 49 (6), 1455-1474. https://doi.org/10.1021/ci900056c.

Doerr, S., Harvey, M.J., Noé, F., Fabritiis, G.D., 2016. HTMD: high-throughput molecular dynamics for molecular discovery. J. Chem. Theory Comput. 12 (4), 1845-1852. https://doi.org/10.1021/acs.jctc.6b00049.

Dunbrack, R.L., Karplus, M., 1994. Conformational analysis of the backbone-dependent rotamer preferences of protein sidechains. Nat. Struct. Mol. Biol. 1 (5), 334-340. https://doi.org/10.1038/nsb0594-334.

Eberhardt, J., Santos-Martins, D., Tillack, A.F., Forli, S., 2021. AutoDock vina 1.2.0: new docking methods, expanded force field, and python bindings. J. Chem. Inf. Model. 61 (8), 3891-3898. https://doi.org/10.1021/acs.jcim.1c00203.

Fisher, R.A., 1992. Statistical Methods for Research Workers. Springer, New York, pp. 66-70. https://doi.org/10.1007/978-1-4612-4380-9_6.

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., Li, X., Caricato, M., Marenich, A.V., Bloino, J., Janesko, B.G., Gomperts, R., Mennucci, B., Hratchian, H. P., Ortiz, J.V., Izmaylov, A.F., Sonnenberg, J.L., Williams-Young, D., Ding, F., Lipparini, F., Egidì, F., Goings, J., Peng, B., Petrone, A., Henderson, T., Ranasinghe, D., Zakrzewski, V.G., Gao, J., Rega, N., Zheng, G., Liang, W., Hada, M., Ehara, M., Toyota, K., Fukuda, R., Hasegawa, J., Ishida, M., Nakajima, T., Honda, Y., Kitao, O., Nakai, H., Vreven, T., Throssell, K., Montgomery Jr., J.A., Peralta, J.E., Ogliaro, F., Bearpark, M.J., Heyd, J.J., Brothers, E.N., Kudin, K.N., Staroverov, V.N., Keith, T.A., Kobayashi, R., Normand, J., Raghavachari, K., Rendell, A.P., Burant, J. C., Iyengar, S.S., Tomasi, J., Cossi, M., Millam, J.M., Klene, M., Adamo, C., Cammi, R., Ochterski, J.W., Martin, R.L., Morokuma, K., Farkas, O., Foresman, J.B., Fox, D.J., 2016. Gaussian Inc ~ 16 Revision C.01. Gaussian Inc, Wallingford CT.

Gill, R.J., Ramos-Rodriguez, O., Raine, N.E., 2012. Combined pesticide exposure severely affects individual- and colony-level traits in bees. Nature 491 (7422), 105-108. https://doi.org/10.1038/nature11585.

Halgren, T.A., 1996. Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94. J. Comput. Chem. 17 (5-6), 490-519. https://doi.org/10.1002/(sici)1096-987x(199604)17:5<6490::aid-jccl>3.0:co-2-p.

Han, W., Tian, Y., Shen, X., 2018. Human exposure to neonicotinoid insecticides and the evaluation of their potential toxicity: an overview. Chemosphere 192, 59-65. https://doi.org/10.1016/j.chemosphere.2017.10.149.

Henry, M., Béguin, M., Requier, F., Rollin, O., Odoux, J.-F., Aupinel, P., Aptel, J., Tchamitchian, S., Decourtye, A., 2012. A common pesticide decreases foraging success and survival in honey bees. Science 336 (6079), 348-350. https://doi.org/10.1126/science.1215039.

Huang, J., Rauscher, S., Nawrocki, G., Ran, T., Feig, M., de Groot, B.L., Grubmüller, H., MacKerell, A.D., 2016. CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nat. Methods 14 (1), 71-73. https://doi.org/10.1038/nmeth.4067.

Huang, S.-Y., Zou, X., 2010. Advances and challenges in protein-ligand docking. Int. J. Mol. Sci. 11 (8), 3016-3034. https://doi.org/10.3390/ijms11083016.

Jeschke, P., Nauen, R., Schindler, M., Elbert, A., 2011. Overview of the status and global strategy for neonicotinoids. J. Agric. Food Chem. 59 (7), 2897-2908. https://doi.org/10.1021/jf101303g.

Kaélé, L., Skeel, R., Bhandarkar, M., Brunner, R., Gursoy, A., Krawetz, N., Phillips, J., Shinozaki, A., Varadarajan, K., Schulten, K., 1999. NAMD2: greater scalability for parallel molecular dynamics. J. Comput. Phys. 151 (1), 283-312. https://doi.org/10.1006/jcph.1999.6201.
Kitchen, D.B., Decornez, H., Furr, J.R., Bajorath, J., 2004. Docking and scoring in virtual screening for drug discovery: methods and applications. Nat. Rev. Drug Discov. 3 (11), 935-949. https://doi.org/10.1038/ndr1549.

Landrum, G., 2021. Rdkit (q1 2021) release. (https://doi.org/10.5281/ZENODO.4750957), (https://zenodo.org/record/4750957).

Mansouri, K., Abdelaziz, A., Rybacka, A., Roncaglioni, A., Tropsha, A., Varnek, A., Zakharov, A., Worth, A., Richard, A.M., Grulke, C.M., Trisciuzzi, D., Fourches, D., Horvath, D., Benfenati, E., Muratov, E., Wedebye, E.B., Grisoni, F., Mangiatordi, G. F., Incisivo, G.M., Hong, H., Ng, H.W., Tetko, I.V., Balabin, I., Kancherla, J., Shen, J., Burton, J., Nicklaus, M., Cassotti, M., Nikolov, N.G., Nicolotti, O., Andersson, P.L., Zang, Q., Politi, R., Beger, R.D., Todeschini, R., Huang, R., Farag, S., Rosenberg, S.A., Slavov, S., Hu, X., Judson, R.S., 2016. CERAPP: collaborative estrogen receptor activity prediction project. Environ. Health Perspect. 124 (7), 1023-1033. https://doi.org/10.1289/ehp.1510267.

Marrs, T.C., 2012. Toxicology of insecticides to mammals. Pest Manag. Sci. 68 (10), 1332-1336. https://doi.org/10.1002/ps.3362.

Morris, G.M., Lim-Wilby, M., 2008. Molecular docking. In: Methods in Molecular Biology. Humana Press, pp. 365-382. https://doi.org/10.1007/978-1-59745-177-2_19.

Noviello, C.M., Gharpure, A., Mukhtasimova, N., Cabuco, R., Baxter, L., Borek, D., Sine, S.M., Hibbs, R.E., 2021. Structure and gating mechanism of the α7 nicotinic acetylcholine receptor. Cell 184 (8), 2121-2134.e13. https://doi.org/10.1016/j.cell.2021.02.049.

O'Boyle, N.M., Banck, M., James, C.A., Morley, C., Vandermeersch, T., Hutchison, G.R., 2011. Open babel: an open chemical toolbox. J. Cheminf. 3 (1) https://doi.org/10.1186/1758-2946-3-33.

Olsson, M.H.M., Søndergaard, C.R., Rostkowski, M., Jensen, J.H., 2011. PROPKA3: consistent treatment of internal and surface residues in empirical pk<sub>a</sub> predictions. J. Chem. Theory Comput. 7 (2), 525-537. https://doi.org/10.1021/ct100578z.

Phillips, J.C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., Chipot, C., Skeel, R.D., Kalé, L., Schulten, K., 2005. Scalable molecular dynamics with NAMD. J. Comput. Chem. 26 (16), 1781-1802. https://doi.org/10.1002/jcc.20289.

Phillips, J.C., Hardy, D.J., Maia, J.D.C., Stone, J.E., Ribeiro, J.V., Bernardi, R.C., Buch, R., Fiorin, G., Hénin, J., Jiang, W., McGreevy, R., Melo, M.C.R., Radak, B.K., Skeel, R.D., Singharoy, A., Wang, Y., Roux, B., Aksimentiev, A., Luthey-Schulten, Z., Kalé, L.V., Schulten, K., Chipot, C., Tajkhorshid, E., 2020. Scalable molecular dynamics on cpu and gpu architectures with namd. J. Chem. Phys. 153 (4), 044130 https://doi.org/10.1063/5.0014475.

Pisa, L.W., Amaral-Rogers, V., Belzunces, L.P., Bonmatin, J.M., Downs, C.A., Goulson, D., Kreutzweiser, D.P., Krupke, C., Liess, M., McField, M., Morrissey, C.A., Noome, D.A., Settele, J., Simon-Delso, N., Stark, J.D., derSluijs, J.P.V., Dyck, H.V., Wiemers, M., 2014. Effects of neonicotinoids and fipronil on non-target invertebrates. Environ. Sci. Pollut. Res. 22 (1), 68-102. https://doi.org/10.1007/s11356-014-3471-x.

Ravindranath, P.A., Forli, S., Goodsell, D.S., Olson, A.J., Sanner, M.F., 2015. Autodockfr: advances in protein-ligand docking with explicitly specified binding site flexibility. PLoS Comput. Biol. 11 (12), e1004586 https://doi.org/10.1371/journal.pcbi.1004586.

Rundlöf, M., Andersson, G.K.S., Bommarco, R., Fries, I., Hederström, V., Herbertsson, L., Jonsson, O., Klatt, B.K., Pedersen, T.R., Yourstone, J., Smith, H.G., 2015. Seed coating with a neonicotinoid insecticide negatively affects wild bees. Nature 521 (7550), 77-80. https://doi.org/10.1038/nature14420.

Schneider, G., Fechner, U., 2005. Computer-based de novo design of drug-like molecules. Nat. Rev. Drug Discov. 4 (8), 649-663. https://doi.org/10.1038/ndr1799.

Schrödinger release, 2021. epik, Schrödinger, LLC, New York, NY.

Sieg, J., Flachsenberg, F., Rarey, M., 2019. In need of bias control: evaluating chemical data for machine learning in structure-based virtual screening. J. Chem. Inf. Model. 59 (3), 947-961. https://doi.org/10.1021/acs.jcim.8b00712.

Spurny, R., Debaveye, S., Farinha, A., Veys, K., Vos, A.M., Gossas, T., Atack, J., Bertrand, S., Bertrand, D., Danielson, U.H., Tresadern, G., Ulens, C., 2015. Molecular blueprint of allosteric binding sites in a homologue of the agonist-binding domain of the α7 nicotinic acetylcholine receptor. Proc. Natl. Acad. Sci. USA 112 (19). https://doi.org/10.1073/pnas.1418289112.

The open babel package, 2020. Version 3.1.0. (http://openbabel.org).

Trott, O., Olson, A.J., 2009. AutoDock vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. https://doi.org/10.1002/jcc.21334 (NA-NA).

Tsvetkov, N., Samson-Robert, O., Sood, K., Patel, H.S., Malena, D.A., Gajiwala, P.H., Maciukiewicz, P., Fournier, V., Zayed, A., 2017. Chronic exposure to neonicotinoids

reduces honey bee health near corn crops. Science 356 (6345), 1395-1397. https://doi.org/10.1126/science.aam7470.

van der Sluijs, J.P., Amaral-Rogers, V., Belzunces, L.P., van Lexmond, M.F.I.J.B., Bonmatin, J.-M., Chagnon, M., Downs, C.A., Furlan, L., Gibbons, D.W., Giorio, C., Girolami, V., Goulson, D., Kreutzweiser, D.P., Krupke, C., Liess, M., Long, E., McField, M., Mineau, P., Mitchell, E.A.D., Morrissey, C.A., Noome, D.A., Pisa, L., Settele, J., Simon-Delso, N., Stark, J.D., Tapparo, A., Dyck, H.V., van Praagh, J., Whitehorn, P.R., Wiemers, M., 2014. Conclusions of the worldwide integrated assessment on the risks of neonicotinoids and fipronil to biodiversity and ecosystem functioning. Environ. Sci. Pollut. Res. 22 (1), 148-154. https://doi.org/10.1007/s11356-014-3229-5.

Vanommeslaeghe, K., Hatcher, E., Acharya, C., Kundu, S., Zhong, S., Shim, J., Darian, E., Guvench, O., Lopes, P., Vorobyov, I., Mackerell, A.D., 2010. CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 31, 671-690. https://doi.org/10.1002/jcc.21367.

Vanommeslaeghe, K., MacKerell, A.D., 2012. Automation of the CHARMM general force field (CGenFF) i: bond perception and atom typing. J. Chem. Inf. Model. 52 (12), 3144-3154. https://doi.org/10.1021/ci300363c.

Vanommeslaeghe, K., Raman, E.P., MacKerell, A.D., 2012. Automation of the CHARMM general force field (CGenFF) II: assignment of bonded parameters and partial atomic charges. J. Chem. Inf. Model. 52 (12), 3155-3168. https://doi.org/10.1021/ci3003649.

Verdonk, M.L., Cole, J.C., Hartshorn, M.J., Murray, C.W., Taylor, R.D., 2003. Improved protein-ligand docking using GOLD. Proteins: Struct. Funct. Bioinf. 52 (4), 609-623. https://doi.org/10.1002/prot.10465.

Walsh, R.M., Roh, S.-H., Gharpure, A., Morales-Perez, C.L., Teng, J., Hibbs, R.E., 2018. Structural principles of distinct assemblies of the human α4/β2 nicotinic receptor. Nature 557 (7704), 261-265. https://doi.org/10.1038/s41586-018-0081-7.

Wang, R., Fang, X., Lu, Y., Yang, C.-Y., Wang, S., 2005. The pdbbind database: methodologies and updates. J. Med. Chem. 48 (12), 4111-4119. https://doi.org/10.1021/jm048957q.

Wang, A., Mahai, G., Wan, Y., Yang, Z., He, Z., Xu, S., Xia, W., 2020. Assessment of imidacloprid related exposure using imidacloprid-olefin and desnitro-imidacloprid: neonicotinoid insecticides in human urine in wuhan, china. Environ. Int. 141, 105785 https://doi.org/10.1016/j.envint.2020.105785.

Wang, Z., Sun, H., Yao, X., Li, D., Xu, L., Li, Y., Tian, S., Hou, T., 2016. Comprehensive evaluation of ten docking programs on a diverse set of protein-ligand complexes: the prediction accuracy of sampling power and scoring power. Phys. Chem. Chem. Phys. 18 (18), 12964-12975. https://doi.org/10.1039/c6cp01555g.

Warren, G.L., Andrews, C.W., Capelli, A.-M., Clarke, B., LaLonde, J., Lambert, M.H., Lindvall, M., Nevins, D., Semus, S.F., Senger, S., Tedesco, G., Wall, I.D., Woolven, J. M., Peishoff, C.E., Head, M.S., 2005. A critical assessment of docking programs and scoring functions. J. Med. Chem. 49 (20), 5912-5931. https://doi.org/10.1021/jm050362n.

Whitehorn, P.R., O'Connor, S., Wackers, F.L., Goulson, D., 2012. Neonicotinoid pesticide reduces bumble bee colony growth and queen production. Science 336 (6079), 351-352. https://doi.org/10.1126/science.1215025.

Woodcock, B.A., Isaac, N.J.B., Bullock, J.M., Roy, D.B., Garthwaite, D.G., Crowe, A., Pywell, R.F., 2016. Impacts of neonicotinoid use on long-term population changes in wild bees in england. Nat. Commun. 7 (1) https://doi.org/10.1038/ncomms12459.

Wu, G., Robertson, D.H., Brooks, C.L., Vieth, M., 2003. Detailed analysis of grid-based molecular docking: a case study of CDOCKER—a CHARMM-based MD docking algorithm. J. Comput. Chem. 24 (13), 1549-1562. https://doi.org/10.1002/jcc.10306.

Xiao, M., 2018. Getbox-pymol-plugin, v20180204. (https://github.com/MengwuXiao/GetBox-PyMOL-Plugin).

Xue, Q., Liu, X., Russell, P., Li, J., Pan, W., Fu, J., Zhang, A., 2022. Evaluation of the binding performance of flavonoids to estrogen receptor alpha by autodock, autodock vina and surflex-dock. Ecotoxicol. Environ. Saf. 233, 113323 https://doi.org/10.1016/j.eecoenv.2022.113323.

Zhang, H., Zhang, R., Zeng, X., Wang, X., Wang, D., Jia, H., Xu, W., Gao, Y., 2022. Exposure to neonicotinoid insecticides and their characteristic metabolites: association with human liver cancer. Environ. Res. 208, 112703 https://doi.org/10.1016/j.envres.2022.112703.

Zhao, H., Gafirsch, A., 2013. Discovery of zap70 inhibitors by high-throughput docking into a conformation of its kinase domain generated by molecular dynamics. Bioorg. Med. Chem. Lett. 23 (20), 5721-5726. https://doi.org/10.1016/j.bmcl.2013.08.009.