Network Toxicology and Molecular Docking to Investigate the Non-AChE Mechanisms of Organophosphate-Induced Neurodevelopmental Toxicity

Organophosphate pesticides (OPs) are toxic substances that contaminate aquatic environments, interfere with the development of the nervous system, and induce Neurodevelopmental Toxicity (NDT) in animals and humans. The canonical mechanism of OP neurotoxicity involves the inhibition of acetylcholinesterase (AChE), but other mechanisms non-AChE are also involved and not fully understood. We used network toxicology and molecular docking to identify molecular targets and toxicity mechanisms common to OPs. Targets related to diazinon-oxon, chlorpyrifos oxon, and paraoxon OPs were predicted using the Swiss Target Prediction and PharmMapper databases. Targets related to NDT were compiled from GeneCards and OMIM databases. In order to construct the protein–protein interaction (PPI) network, the common targets between OPs and NDT were imported into the STRING. Network topological analyses identified EGFR, MET, HSP90AA1, and SRC as hub nodes common to the three OPs. Using the Reactome pathway and gene ontology, we found that signal transduction, axon guidance, cellular responses to stress, and glutamatergic signaling activation play key roles in OP-induced NDT.


Introduction
The widespread use of organophosphate pesticides (OPs) in agriculture still raises concerns about their ecological impact and human health implication [1]. Surface waters and sediments can be contaminated by OPs from rainfall, runoff, improper disposal, or leaching from groundwater. This contamination can affect a wide range of aquatic organisms, such as fish, invertebrates, and algae, as well as non-target terrestrial organisms [2].
There are numerous toxic effects associated with OPs in fish, including hematological disturbances, gill, kidney, and liver alterations [3], oxidative stress, immune disorders, alteration of the intestinal microbiota [4], and behavioral disorders including a reduced predator escape response [5]. In humans, exposure to OPs during pregnancy and postnatal periods increases the risk of autism spectrum disorder [6], impaired IQ, and verbal comprehension [7].
Exposure to OPs poses a high risk of neurotoxicity, especially for the developing nervous system. The effects of this type of exposure can include cognitive, motor, and memory impairments, as well as changes in brain morphology, referred to as Neurodevelopmental Toxicity (NDT) [8].

Chemical Structures of OPs
The 2D structures of the oxon metabolites of the OP pesticides diazinon (DZ), CP, and parathion were retrieved from the PubChem database (h ps://pubchem.ncbi.nlm.nih.gov) accessed on 28 February 2023. We chose the oxon metabolites over of the original compounds due to their direct association with NDT [20,21]. The images of the chemical structures of the organophosphates were generated using Marvin JS [22] and are shown in Figure 1.

Acquisition of OP Targets
The 2D chemical structures of paraoxon (PO), chlorpyrifos oxon (CPO), and diazinon-oxon (DO) were submi ed to the Swiss Target Prediction database (h p://www.swisstargetprediction.ch) accessed on 28 February 2023 and PharmMapper (h p://www.lilab-ecust.cn/pharmmapper/) accessed on 28 February 2023 to predict target genes. We then entered the predicted target genes into the UniProt database (h ps://www.uniprot.org/) accessed on 28 February 2023 to retrieve the standard gene names, and duplicate entries were eliminated.

Prediction of NDT Targets
The GeneCards database (h ps://www.genecards.org/) accessed on 28 February 2023 and Online Mendelian Inheritance in Man (OMIM, h p://omim.org/) accessed on 28 February 2023 were used to predict potential NDT targets. In February 2023, the related targets were collected using the keywords "neurodevelopmental abnormalities", "neurotoxicity", and "neurodevelopmental disorders", along with the Homo sapiens species. The UniProt database was utilized to retrieve the standard gene names, and after merging the predicted targets for the three keywords from both databases, duplicate entries were removed.

Venn Analysis and Construction of Protein-Protein Interaction (PPI) Network
In order to identify the intersection genes as potential targets of each OP related to NDT, a Venn diagram was created by using Draw Venn Diagram tool (h p://bioinformatics.psb.ugent.be/webtools/Venn/) accessed on 28 February 2023.
Protein-protein interaction (PPI) networks were constructed using the STRING database (h ps://string-db.org/) accessed on 28 February 2023. A minimum score of 0.7 was set, indicating a high level of confidence, and only experimental data were selected as the source of active interaction. The PPI network analyzes direct interactions between the inpu ed proteins and generates a network of interactions where proteins are represented as nodes and the interactions between them are represented as edges.

Acquisition of OP Targets
The 2D chemical structures of paraoxon (PO), chlorpyrifos oxon (CPO), and diazinon-oxon (DO) were submitted to the Swiss Target Prediction database (http://www.swisstargetprediction. ch) accessed on 28 February 2023 and PharmMapper (http://www.lilab-ecust.cn/pharmmapper/) accessed on 28 February 2023 to predict target genes. We then entered the predicted target genes into the UniProt database (https://www.uniprot.org/) accessed on 28 February 2023 to retrieve the standard gene names, and duplicate entries were eliminated.

Prediction of NDT Targets
The GeneCards database (https://www.genecards.org/) accessed on 28 February 2023 and Online Mendelian Inheritance in Man (OMIM, http://omim.org/) accessed on 28 February 2023 were used to predict potential NDT targets. In February 2023, the related targets were collected using the keywords "neurodevelopmental abnormalities", "neurotoxicity", and "neurodevelopmental disorders", along with the Homo sapiens species. The UniProt database was utilized to retrieve the standard gene names, and after merging the predicted targets for the three keywords from both databases, duplicate entries were removed.

Venn Analysis and Construction of Protein-Protein Interaction (PPI) Network
In order to identify the intersection genes as potential targets of each OP related to NDT, a Venn diagram was created by using Draw Venn Diagram tool (http://bioinformatics.psb. ugent.be/webtools/Venn/) accessed on 28 February 2023.
Protein-protein interaction (PPI) networks were constructed using the STRING database (https://string-db.org/) accessed on 28 February 2023. A minimum score of 0.7 was set, indicating a high level of confidence, and only experimental data were selected as the source of active interaction. The PPI network analyzes direct interactions between the inputted proteins and generates a network of interactions where proteins are represented as nodes and the interactions between them are represented as edges.

Topological Analysis of PPI Networks
The PPI networks were imported into Cytoscape 3.9.1. The cytoHubba plugin was utilized to compute the top 10 hub nodes using the Maximal Clique Centrality (MCC) method. The Analyze Network tool was used for topological analyses of the network, including degree (k), betweenness centrality (BC), clustering coefficient, closeness centrality (CC), and average shortest pathway (ASPL).
The measure k represents the number of connections a node has with other nodes in the network, i.e., it quantifies the number of edges that connect to it. Nodes with a high k value are important for information transmission in the network as they have more control over the flow of information. On the other hand, BC is a measure that indicates a node's ability to connect other nodes in the network. A high BC value indicates that the node is important for communication between different parts of the network. In turn, the CC measure indicates how many direct neighbors of a node are interconnected, and a high CC value for a node indicates its involvement in a strongly connected cluster. Lastly, the ASPL measure indicates the average shortest path length between a node and all other nodes in the network [23].

Enrichment Analysis
An enrichment analysis was performed using the Consensus PathDB database (http: //cpdb.molgen.mpg.de/) accessed on 28 February 2023 to explore the Gene Ontology (GO) and Reactome Pathway. The Biological Process (BP), Molecular Function (MF), and Cellular Component (CC) categories of GO were chosen. GO terms and Reactome Pathways with a p-value < 0.01 were considered statistically significant. The top 15 GO terms for BP, MF, and CC, as well as the top 30 Reactome pathways, were selected for each OP. Additionally, terms common to all three OPs that were directly related to the Nervous System superpathway, as classified by Reactome.org, were examined. GO and Reactome pathway graphs were generated with the ggplot2 package in RStudio.

Molecular Docking
A molecular docking analysis was conducted to explore the potential interactions between OPs and their hub nodes. The analysis followed the following steps: Docking simulation was performed using Autodock Vina [24], using a grid box with dimensions of 25 × 25 × 25 Å. The xyz coordinates obtained from Discovery Studio were used to select the center of mass for each macromolecule ( Table 1). The poses were selected based on the lowest root-mean-square deviation (RMSD) values, with a maximum threshold of 2.0 Å. (d) The analysis of the binding and generation of images was carried out using Discovery Studio 2021. The distance criterion between ligands and amino acid residues was established as <3.3 Å for hydrogen bonds [25] and <6.0 Å for π-π, π-alkyl and π-sigma interactions [26].

Candidate Targets for OP-Induced NDT
After removing the duplicates, 388 targets were found for DZO, 343 for CPO, and 387 for PO in the Swiss Target Prediction and PharmMapper databases. Additionally, 3060 NDT-related targets were obtained from the Gene Cards and OMIM databases after removing the duplicates. Using a Venn diagram to identify the intersections of targets for each OP and NDT, we found 187, 121, and 181 targets for DZO, CPO, and PO, respectively. These targets were used to construct a PPI network using the STRING database, and the data were imported into Cytoscape for visualization and analysis of the hub nodes.
The analysis of PPI networks revealed 81 nodes and 113 edges for the DZO-NDT network, 37 nodes and 50 edges for the CPO-NDT network, and 71 nodes and 103 edges for the PO-NDT network ( Figure 2). The top 10 hub nodes for each network were identified using the MCC algorithm in the cytoHubba plugin of Cytoscape. Topological measurements were obtained using the Network Analyzer tool. HSP90AA1, SRC, MET, and EGFR were hub nodes that were common among the three OPs. All hub nodes are listed in Table 2, along with their topological measurements. The measurements for the entire network are available in the supplementary materials (Tables S1-S3).      Figure 3).
We observed that the top 15 GO terms for BP, MF, and CC were similar for DZO, CPO, and PO. The results indicated that the BP terms were mainly related to cellular response to stress, which could be a result of chemical exposure, as well as cellular processes such as communication, signal transduction, and cell death, along with processes involved in the development of an anatomical structure or organism.
Regarding the MF terms, we found that most targets for the three OPs were associated with binding, which refers to the interaction between one molecule and specific sites of another molecule. The results of the CC terms suggested most targets were found within the intracellular milieu, organelles, plasma membrane, and neurons.  (Figure 3).
We observed that the top 15 GO terms for BP, MF, and CC were similar for DZO, CPO, and PO. The results indicated that the BP terms were mainly related to cellular response to stress, which could be a result of chemical exposure, as well as cellular processes such as communication, signal transduction, and cell death, along with processes involved in the development of an anatomical structure or organism.
Regarding the MF terms, we found that most targets for the three OPs were associated with binding, which refers to the interaction between one molecule and specific sites of another molecule. The results of the CC terms suggested most targets were found within the intracellular milieu, organelles, plasma membrane, and neurons. (a)

Reactome Pathway Analysis
Based on Reactome pathway analysis, 419 pathways were identified for DZO, 254 for CPO, and 358 for PO (p-value < 0.01). Figure 4 shows the top 30 Reactome pathways for each OP based on the p-value.

Reactome Pathway Analysis
Based on Reactome pathway analysis, 419 pathways were identified for DZO, 254 for CPO, and 358 for PO (p-value < 0.01). Figure 4 shows the top 30 Reactome pathways for each OP based on the p-value.  Reactome pathway is a database that encompasses 2546 pathways organized into 28 superpathways (e.g., signal transduction, disease, immune system, and developmental biology) and their subdivisions [27]. According to the analyses, signal transduction was identified as the most important superpathway for DZO, CPO, and PO, which included pathways related to signaling by receptor tyrosine kinases, MAPK family signaling cascades, and intracellular signaling by second messengers. Additionally, all three OPs Reactome pathway is a database that encompasses 2546 pathways organized into 28 superpathways (e.g., signal transduction, disease, immune system, and developmental biology) and their subdivisions [27]. According to the analyses, signal transduction was identified as the most important superpathway for DZO, CPO, and PO, which included pathways related to signaling by receptor tyrosine kinases, MAPK family signaling cascades, and intracellular signaling by second messengers. Additionally, all three OPs shared the axon guidance pathway, which is directly related to nervous system development, suggesting that OPs interfere with neurogenesis.
We also evaluated the Reactome pathways for DZO, CPO and PO that showed significant values (p-value < 0.01) and were classified under the Neuronal system superpathway to identify the key processes in the nervous system affected by OP exposure. Activation of NMDA receptors and postsynaptic events was the common pathway for all three OPs, indicating the activation of N-methyl-D-aspartate receptors (NMDAR), an ionotropic glutamate receptor.

Molecular Docking
We used molecular docking to predict the interactions between each OP and their respective hub nodes (presented in Table 2). A lower binding energy between the compound and the target indicates a higher affinity. The binding energies between the OPs and macromolecules ranged from −7.2 to −4.2 kcal/mol ( Figure 5). Among the targets shared by the OPs, HSP90AA1 had the lowest binding energy with the OPs, followed by MET, EGFR, and SRC. shared the axon guidance pathway, which is directly related to nervous system development, suggesting that OPs interfere with neurogenesis. We also evaluated the Reactome pathways for DZO, CPO and PO that showed significant values (p-value < 0.01) and were classified under the Neuronal system superpathway to identify the key processes in the nervous system affected by OP exposure. Activation of NMDA receptors and postsynaptic events was the common pathway for all three OPs, indicating the activation of N-methyl-D-aspartate receptors (NMDAR), an ionotropic glutamate receptor.

Molecular Docking
We used molecular docking to predict the interactions between each OP and their respective hub nodes (presented in Table 2). A lower binding energy between the compound and the target indicates a higher affinity. The binding energies between the OPs and macromolecules ranged from −7.2 to −4.2 kcal/mol ( Figure 5). Among the targets shared by the OPs, HSP90AA1 had the lowest binding energy with the OPs, followed by MET, EGFR, and SRC. The interactions between the OPs and the common hub nodes, EGFR, HSP90AA1, MET, and SRC, are shown in Figure 6. The OPs interacted with the ATP-binding site within the N-terminal domain of HSP90AA1. The ligands established π-sigma, π-alkyl and π-π interactions with specific amino acid residues. Similarly to PU3, an HSP90 inhibitor, the three OPs interacted with LEU 107 and PHE 138 , and PO also formed hydrogen bonds with TRP 162 [28]. The interactions between OPs and EGFR occurred in the ATP-binding cleft, positioned between the amino-terminal and carboxi-terminal lobes-an allosteric site. These interactions involved the formation of π-sigma, π-alkyl and hydrogen bonds with the amino acid residues. Similar to erlotinib, an EGFR inhibitor, the OPs interacted with MET 769 and THR 766 residues, as well as LYS 721 , recognized as one of the key residues for EGFR biological activity. Furthermore, all three OPs exhibited a ractive forces with ASP 813 residue of the catalytic loop [29,30].
OPs established multiple π-alkyl interactions with MET, in addition to forming πsigma and π-π bonds. Interactions with residues at the ATP binding site were observed, including PHE 1089 and VAL 1092 in the N-lobe, and LEU 1157 in the hinge region of the binding pocket. This region engages with non-competitive ATP inhibitors like Tivantinib [31,32]. The interactions between the OPs and the common hub nodes, EGFR, HSP90AA1, MET, and SRC, are shown in Figure 6. The OPs interacted with the ATP-binding site within the N-terminal domain of HSP90AA1. The ligands established π-sigma, π-alkyl and π-π interactions with specific amino acid residues. Similarly to PU3, an HSP90 inhibitor, the three OPs interacted with LEU 107 and PHE 138 , and PO also formed hydrogen bonds with TRP 162 [28]. The interactions between OPs and EGFR occurred in the ATP-binding cleft, positioned between the amino-terminal and carboxi-terminal lobes-an allosteric site. These interactions involved the formation of π-sigma, π-alkyl and hydrogen bonds with the amino acid residues. Similar to erlotinib, an EGFR inhibitor, the OPs interacted with MET 769 and THR 766 residues, as well as LYS 721 , recognized as one of the key residues for EGFR biological activity. Furthermore, all three OPs exhibited attractive forces with ASP 813 residue of the catalytic loop [29,30].
The interactions between the OPs and the SRC SH2 domain occurred through hydrogen bonding, π-sigma, and π-alkyl interactions. DZO and CPO established hydrogen bonds with the LYS 62 residue. Furthermore, DZO formed a hydrogen bond with ARG 14 at the phosphotyrosine binding site. PO engaged in hydrogen bonding with the Leu 96 , Gly 95 , and Tyr 89 residues within the specificity pocket [33].
Detailed results of the molecular docking analysis can be found in the Supplementary Materials (Tables S4 and S5 and Figure S1).

Discussion
The classic mechanism of neurotoxicity of organophosphates is the inhibition of the AChE activity, leading to the accumulation of the neurotransmi er acetylcholine in the synaptic cleft and consequently, overstimulation of cholinergic receptors. However, numerous studies have provided evidence that multiple mechanisms can contribute to the neurotoxicity of organophosphates, particularly during neurodevelopment [12].
Reactome pathway and Gene ontology enrichment analyzes were performed to identify potential mechanisms in the NTD of OP. The results indicate that these compounds impair signal transduction and axon guidance, disrupt cellular responses to stress, and activate NMDAR. In the OP-induced NTD PPI network, HSP90AA1, EGFR, MET, and SRC were identified as common hub nodes for DZO, CPO, and PO. For be er understanding, the discussion will be divided into several topics according to the mechanism of toxicity.

Signal Transduction
Signal transduction is the process by which extracellular messengers bind to transmembrane receptors and provide information that triggers cellular responses, such as biochemical and biological transformations, as well as gene expression alterations [34]. According to Reactome pathway enrichment analyses, the major signal transduction pathways affected by organophosphates include the receptor Tyrosine Kinase (RTK) pathway, the mitogen-activated protein kinase (MAPK) pathway and signaling by second messengers. It is important to note that these pathways are interconnected, since RTK activation triggers phosphorylation and the recruitment of effector proteins, initiating signaling cascades involving MAPKs, PI3K/Akt, phospholipase C-PKC, and STAT [35]. OPs established multiple π-alkyl interactions with MET, in addition to forming πsigma and π-π bonds. Interactions with residues at the ATP binding site were observed, including PHE 1089 and VAL 1092 in the N-lobe, and LEU 1157 in the hinge region of the binding pocket. This region engages with non-competitive ATP inhibitors like Tivantinib [31,32]. The interactions between the OPs and the SRC SH2 domain occurred through hydrogen bonding, π-sigma, and π-alkyl interactions. DZO and CPO established hydrogen bonds with the LYS 62 residue. Furthermore, DZO formed a hydrogen bond with ARG 14 at the phosphotyrosine binding site. PO engaged in hydrogen bonding with the Leu 96 , Gly 95 , and Tyr 89 residues within the specificity pocket [33].
Detailed results of the molecular docking analysis can be found in the Supplementary Materials (Tables S4 and S5 and Figure S1).

Discussion
The classic mechanism of neurotoxicity of organophosphates is the inhibition of the AChE activity, leading to the accumulation of the neurotransmitter acetylcholine in the synaptic cleft and consequently, overstimulation of cholinergic receptors. However, numerous studies have provided evidence that multiple mechanisms can contribute to the neurotoxicity of organophosphates, particularly during neurodevelopment [12].
Reactome pathway and Gene ontology enrichment analyzes were performed to identify potential mechanisms in the NTD of OP. The results indicate that these compounds impair signal transduction and axon guidance, disrupt cellular responses to stress, and activate NMDAR. In the OP-induced NTD PPI network, HSP90AA1, EGFR, MET, and SRC were identified as common hub nodes for DZO, CPO, and PO. For better understanding, the discussion will be divided into several topics according to the mechanism of toxicity.

Signal Transduction
Signal transduction is the process by which extracellular messengers bind to transmembrane receptors and provide information that triggers cellular responses, such as biochemical and biological transformations, as well as gene expression alterations [34]. According to Reactome pathway enrichment analyses, the major signal transduction pathways affected by organophosphates include the receptor Tyrosine Kinase (RTK) pathway, the mitogen-activated protein kinase (MAPK) pathway and signaling by second messengers. It is important to note that these pathways are interconnected, since RTK activation triggers phosphorylation and the recruitment of effector proteins, initiating signaling cascades involving MAPKs, PI3K/Akt, phospholipase C-PKC, and STAT [35].
In vitro studies have explored how OP exposure affects signaling pathways. For instance, exposure of cell cultures to CP and monocrotophos activated the MAPK ERK1/2, JNK, and p38 pathways, thereby inducing apoptosis. Activation of these pathways were associated with reactive oxygen species (ROS) generation and oxidative stress, which may contribute to mitochondrial damage and cell death [36][37][38]. It is plausible that these pathways may also be influenced by the interaction of OPs with RTKs.
Based on our literature review, no experimental data were found, indicating that OPs cause neurodevelopmental abnormalities via EGFR, MET, or non-receptor tyrosine kinase protein SRC. However, CP has previously been shown to increase migration and invasion of breast cancer cells through c-SRC pathway activation, amplifying downstream AKT and p-38 signaling [39]. Additionally, CP has promoted the growth of human colorectal adenocarcinoma H508 cells through increased EGFR/ERK1/2 signaling [40]. Although activation of the EGFR/ERK1/2 pathway generally stimulates cell growth, it can also promote neuronal cell death [41,42].

Axon Guidance
Axon guidance is the process by which neurons send their axons toward their correct targets for synaptic formation, and it is crucial during neural circuit development. Growing axons have a structure at their tips known as the growth cone, which contains receptors that respond to extracellular messengers that attract or repel the axon, guiding them to their destination [43]. The hub nodes EGFR, MET, and SRC play important roles in axon guidance.
Axonal growth is a fundamental process in the formation of neural circuits and the development of neuronal connectivity. Previous studies have reported impairment of axonal growth by OPs through noncholinergic mechanisms. For example, low concentrations of CP and CPO inhibited axonal outgrowth in cultures of rat sympathetic neurons [50]. Moreover, exposure to CPO during neurodevelopment disrupted axonal outgrowth in zebrafish [51]. These results showed that CP and CPO did not cause axon retraction or a decrease in the number of axons, but rather reduced growth rate compared to the untreated control group. Furthermore, no morphological changes in the axonal growth cone were observed. Experimental data suggest that OPs altered axonal growth signaling [50], similar to the findings in our study that OPs impair axon guidance.

Cellular Response to Stress
Exposure to environmental contaminants during prenatal development can trigger cellular response to stress that is closely related to heat shock and oxidative stress systems [52]. The HSP90AA1 gene encodes the Heat shock protein (Hsps) HSP90, which acts as a chaperone to refold denatured proteins in the presence of stressors such as heat, oxidative stress, hypoxia, and exposure to cytotoxic agents. In the absence of stress, chaperones facilitate protein refolding, stabilization, translocation, and degradation, preventing protein aggregation. In neurodevelopment, Hsps regulate pathways related to cell growth and migration, such as the PI3K/Akt signaling pathway, and they are crucial mediators for axon guidance [53].
Cellular stress increases the expression of HSP90 and diverts its functions to cope with stress responses. This functional diversion or inhibition of HSP90 can alter the activity of various proteins and result in increased mutations, leading to neurodevelopmental disorders [54]. The increased expression of HSP90 has been observed in the liver, muscle, kidneys, and spleen of fish exposed to CP, indicating that OPs exposure can induce cellular response to stress [55,56]. Although these data were not obtained from nervous system cells, common carp exposed to CP showed increased expression of HSP90 and hypoactivity, which may indicate neurotoxicity [55].
There are several ways cells respond to stress, including activating of survival pathways or apoptosis [57]. There is still much to be clarified about the pathways by which OPs induce apoptosis, but they are known to occur during vulnerable windows of neurodevelopment and is likely not related to AChE inhibition [58]. However, it has been observed that the AKT pathway and oxidative stress are mechanisms by which OPs induce neural apoptosis, but independently of each other [59]. Therefore, the relationship between OPs and HSP90 may be a pathway to understanding how these processes occur.

Activation of NMDAR
The ionotropic glutamate receptors N-Methyl-D-Aspartate receptors (NMDAR) are calcium-permeable channels. The synaptic activity of NMDAR contributes to axonal and dendritic growth, as well as the maturation of glutamatergic synapses [60]. These processes are crucial for the development of cognitive functions, neuronal plasticity, memory, and learning [61]. Excessive activation of NMDAR can increase calcium influx, leading to cellular damage and neuronal death, likely mediated by the PKC/ERK pathway [62,63].
During exposure to high concentrations of OPs, ionotropic AMPA and NMDA glutamate receptors are activated due to increased excitatory signaling following AChE inhibition and excessive stimulation of muscarinic receptors, leading to seizures and neuronal death [64]. However, the mechanisms related to low concentrations of these chemicals have not been fully elucidated yet.
Studies indicate that glutamate-mediated excitotoxicity seems to be more related to the toxicity effects of CP than DZ. An NMDAR antagonist, MK-801, attenuated CP toxicity, but not DZ toxicity, in cortical neuronal and glial cultures [11]. Furthermore, in neonatal rats exposed to CP and DZ, both substances upregulated the expression of NMDAR subunits, but CP induced more significant changes than DZ [65].
In our study, we did not assess the direct activation of NMDAR by OPs, nor the genes encoding the subunits of this receptor. However, the activity of NMDAR can be regulated by SRC [66], EGFR [67], and MET [68], which have been identified as hub nodes in the PPI network of OPs. Furthermore, studies with chronic exposure to low doses of OPs is a more realistic way to assess environmental exposure, and further studies are needed to evaluate the relationship between this exposure and NMDAR activity, considering other biomarkers beyond the genes encoding the subunits of this receptor.

Conclusions
Using network toxicology and molecular docking analyses, this study identified the hub nodes HSP90AA1, EGFR, MET, and SRC as potential targets for OPs. Additionally, sig-nal transduction, axon guidance, cellular response to stress, and the activation of NMDAR were found to be key pathways involved in OP-induced NDT. These findings are important for understanding the mechanisms of neurotoxicity of these substances at environmentally relevant concentrations that do not involve cholinergic pathways.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/toxics11080710/s1, Table S1: Topological measurements of nodes in the PPI network of diazinon oxon-induced developmental neurotoxicity, Table S2: Topological measurements of nodes in the PPI network of chlorpyrifos oxon-induced developmental neurotoxicity, Table S3: Topological measurements of nodes in the PPI network of paraoxon-induced developmental neurotoxicity, Table S4: Binding energies of molecular docking between organophosphates and their respective hub nodes, Table S5: Intermolecular interactions of complexes between HSP90AA1, EGFR, MET, and SRC and the organophosphates, Figure S1: Protein-ligand interactions between organophosphates and hub nodes generated using BIOVIA Discovery Studio Visualizer.

Data Availability Statement:
The authors confirm that all data underlying the findings are fully available without restriction. All relevant data are within the paper.