Prediction of cholinergic compounds by machine-learning

Journal Pre-proofs Prediction of cholinergic compounds by machine-learning Sanjeeva J. Wijeyesakere, Daniel M. Wilson, M. Sue Marty PII: DOI: Reference:

S2468-1113(20)30001-3 https://doi.org/10.1016/j.comtox.2020.100119 COMTOX 100119

To appear in:

Computational Toxicology

Received Date: Accepted Date:

6 January 2020 10 January 2020

Please cite this article as: S.J. Wijeyesakere, D.M. Wilson, M.S. Marty, Prediction of cholinergic compounds by machine-learning, Computational Toxicology (2020), doi: https://doi.org/10.1016/j.comtox.2020.100119

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier B.V.

Prediction of cholinergic compounds by machine-learning

Sanjeeva J. Wijeyesakere 1,2, Daniel M. Wilson 1 and M. Sue Marty 1 1

The Dow Chemical Company, Midland MI 48674 2

Currently at FMC Corporation, Newark DE

Correspondence should be addressed to: Daniel Wilson, The Dow Chemical Company, Midland, MI 48674. Tel: (989) 636-0712; E-mail: [email protected]

RUNNING HEAD: Computer models to identify cholinergics

Keywords: Cholinergic system, neurobiology, computational biology, machine learning.

Abstract The cholinergic nervous system plays a central role in biology and medicine and is comprised of the nicotinic (nAChR) and muscarinic receptors (mAChR) and acetylcholinesterase (AChE). Using ~20,000 compounds compiled from publicly available data, we developed random-forest machine-learning models that distinguish inhibitors of the cholinergic system from non-inhibitors by analyzing 2-D structural scaffolds and fingerprints. We also developed parallel models that incorporate binding energetics data from protein-ligand docking. The scaffold/fingerprint-based models were highly sensitive (98.4-99.3% AChRs; 100% AChE) and accurate (96.1-96.8% AChRs; 100% AChE). By further incorporating the binding energies from protein docking, the models identified compounds that interact with either the nAChR or mAChR with 100% sensitivity while retaining high accuracy (≥ 91%). Uniquely, these models can predict whether a compound will or will not inhibit the cholinergic system with 91.1-100% balanced accuracy, respectively. In evaluating robustness, ≤12% of available data were needed for highly sensitive models, suggesting that the models will not change for the foreseeable future. Furthermore, we show that training sets for machine learning models requires sufficient diversity of active and inactive (control) compounds for the resulting models to accurately detect a wide array of active substances. Finally, as an example application, when used to query a large in vitro neuro-signaling database, these models flagged all previously identified cholinergic compounds as well as many additional substances with previously uncharacterized mechanisms. These approaches can be readily adopted to model other neuronal and biological protein targets.

2

Introduction The cholinergic system mediates control of the voluntary, sympathetic and parasympathetic branches of the nervous system and represents a major neuronal signaling pathway (Tiwari et al., 2013; Legay and Mei, 2017). The cholinergic system is well-characterized with three discrete targets (Tiwari et al., 2013; Legay and Mei, 2017). Neurotransmission is mediated by acetylcholine (ACh) (Fig. 1), which interacts with either the nicotinic receptor (nAChR; a cys-loop receptor comprising a ligand-gated sodium channel; Tiwari et al., 2013; Fasoli and Gotti, 2015) or the muscarinic receptor (mAChR, a G-protein coupled receptor; Tiwari et al., 2013; Kruse et al., 2014) within the synapse. These receptors are named for their high-affinity alkaloid ligands, nicotine and muscarine (Bikadi and Simonyi, 2003) (Fig. 1). The enzyme acetylcholinesterase (AChE), inactivates ACh via hydrolysis following its interaction with ACh receptors (Silman and Sussman, 2017). The cholinergic system is of keen interest as a pharmaceutical target (e.g. therapies for Alzheimer’s disease; Verma et al., 2018)), a potential off-target of toxicity for compounds (e.g. tri-o-cresyl phosphate, the causative agent of the Ginger Jake epidemic; Richardson et al., 2013)) and to derive counter-measures for chemical agents (e.g. the AChE inhibitors sarin (Fig 1), soman, etc.; Szinicz, 2005)). Although the cholinergic system is of fundamental importance in biology and medicine, there are currently no computational models that distinguish compounds that inhibit the cholinergic system from those that do not. To address this computational gap, we developed semi-automated algorithms using the Konstanz information miner (KNIME) for big data mining of inhibition data for the cholinergic targets from public databases comprised primarily of in vitro high throughput screening (HTS) data. These databases included government-funded efforts like ToxCast and Tox21 (Dix et al., 2007), which have undertaken in vitro HTS assessments of interactions between large compound libraries and biological targets, as well as ChEMBL. We hypothesized that compounds that inhibit the cholinergic system will share a limited set of structural motifs that distinguish them from inactive control compounds. We characterized the 2-D structural features

3

of our database by screening the active (positive) compounds and inactive (control) compounds using 2-D fingerprinting and scaffolding. Fingerprinting condenses predefined libraries of 2D structures to 1D bitvectors that represent the presence or absence of a chemical feature or functional group. Scaffolding derives the sets of substructures in common within a dataset. The actives versus controls were further characterized by developing protein docking models for each cholinergic target that predict binding energies at the ligand binding site. Finally, we developed a random-forest machine-learning model that assimilates all of the information above and predicts whether a novel compound will or will not inhibit the cholinergic system. To the best of our knowledge, the overall computational approach described herein is unique.

Methods Data compilation To determine the cholinergic potential of novel/uncharacterized compounds, we developed an inhouse profiler using the open-source KNIME (ver. 3.6) development environment (Berthold et al., 2008) as highlighted in Fig. 2a. A database of most known cholinergic compounds (corresponding to 20,536 unique compounds that do or do not interact with the nAChR, mAChR and AChE) was compiled from multiple public databases [ChEMBL (Gaulton et al., 2012; Bento et al., 2014), BindingDB (Chen et al., 2001; Gilson et al., 2016), ZINC (Sterling and Irwin, 2015) and ToxCast (Dix et al., 2007)] along with the published literature. For each compound, its name, identifier (e.g., ChEMBL ID, CAS, ZINC ID, etc.), target receptor and structure (encoded as a SMILES pattern) were recorded together with the activity of the compound against the identified target (encoded as a binary value with ‘Yes’ depicting compounds that bound/inhibited the target protein and ‘No’ representing compounds that did not interact with the target). The compiled data were assessed for accuracy and activity (ability to bind/inhibit the target of interest) via manual inspection of the primary data (e.g., dose-response binding curves from ToxCast) or concordance with the published activity measurements. Additionally, a custom set of scaffolds for surfactants (implemented as a set of SMARTS filters) was used to exclude surfactants from our database of cholinergic compounds. Surfactants were excluded from the models since they represent a class of 4

compounds that could artifactually confound HTS assays via non-specific interactions with the receptor (e.g., denaturation or alterations in receptor conformation/conformational dynamics; Tan and Ting, 2000). For organic compounds present in their salt form, we assumed that the activity would be driven by the larger organic component. Additionally, since the primary focus of our research efforts involve organic compounds, we excluded inorganic compounds as well as those containing heavy metals and atoms (As, B, Sn, Si and Zn) whose 3D-structures could not be reliably parametrized using the AMBER14 force field as implemented within the YASARA (ver. 17.8.15) molecular modeling environment. An analogous database of 1,775 compounds that are known to not interact with the cholinergic system was curated from ToxCast as follows: The list of inactive compounds for each target protein was downloaded along with the associated structures (stored as SMILES strings). These compounds were imported into KNIME and compared with the compiled “active” database. A few substances had mixed results reported (i.e., positive and negative) and therefore, were present as both “active “and “inactive” in our database; these substances were excluded from the “inactive” database (favoring model sensitivity). This process yielded a list of compounds with experimental data indicating that they did not interact with the cholinergic targets. The compiled datasets of inhibitors and control compounds are available in Supplemental Tables S1 and S2, respectively.

KNIME workflows to scaffold compounds We implemented a KNIME workflow to structurally scaffold compounds that do or do not interact with the cholinergic system. This was implemented separately for each database (active cholinergic compounds or inactive controls) as follows: The master cholinergic database was divided into three smaller datasets covering substances that do or do not interact with each of the three major proteins that comprise the cholinergic system: the nAChR (4,481 compounds), mAChR (8,931 compounds) or AChE (7,400 compounds). Each compound was fingerprinted using the MACCS algorithm (Durant et al., 2002), yielding a 166 bit linear bit-vector that described its structural features. Following fingerprinting, the compounds were clustered based on their Tanimoto similarity via agglomerative hierarchical clustering such that 5

compounds with Tanimoto similarity values >0.85 (calculated from the MACCS fingerprint) were included within a cluster. Following assignment of the compounds to a cluster, the maximal common substructure (MCS) for each cluster was calculated using the RDKit MCS algorithm in KNIME, yielding an initial set of structural motifs (scaffolds). These preliminary scaffolds were further refined via a second round of clustering based on their Tanimoto similarity using an agglomerative hierarchical clustering algorithm such that scaffolds sharing Tanimoto similarity values >0.75 were included within a cluster. Each cluster of scaffolds was reanalyzed using the RDKit MCS algorithm to calculate their maximal common substructures, yielding a final set of scaffolds that explains the structural diversity of compounds that do or do not target the cholinergic system. Similar to substances identified as both ‘active’ and ‘inactive’ in the compiled datasets, any scaffolds that were present in both datasets (actives and inactive controls) were deleted from the inactive control dataset but retained within the scaffolds for the active cholinergics to ensure the models were optimized for sensitivity. The final scaffold model yielded sets of unique scaffolds that do or do not target the cholinergic system. Two broad assumptions were made in developing the scaffolding workflow: 1) Scaffolds had to contain more than 4 bonds; and 2) scaffolds comprising a sole benzene ring (in the absence of any other structural motifs) were excluded due to its non-descript nature. An overview of the KNIME workflow used to generate the scaffolds is shown in Figure 2.

Prediction models to assess the cholinergic potential of compounds Assessment of the predictive power of our scaffolding analyses was undertaken within KNIME. A hybrid approach was developed using the composite outcome of two sets of random-forest machinelearning models as covariates: 1) a scaffold-based model that used similarity of a test substance to scaffolds calculated from the active and inactive control compounds; and 2) using the MACCS structural fingerprints of actives and inactives as the sole covariate. When constructing the models, a random subset (10%) of the master database of active and inactive control compounds was withheld as a test set, with the remaining 6

compounds (90% of the dataset) used to train the model. Each model was built 5 times with a different random subset of actives and inactive controls with 10% of the data being withheld to test the model. The mean concordance statistics are reported. This modeling approach is described in greater detail below. For the development of the scaffold-based prediction models, the training set of active and inactive control (non-cholinergic) compounds (90% of the data) were scaffolded via algorithm no. 1 described above (active and inactive control datasets were scaffolded separately). The cholinergic and non-cholinergic (inactive control) scaffolds were then fingerprinted using the MACCS algorithm and matched to the MACCS fingerprints of the withheld active compounds or the database of inactive (non-cholinergic) compounds within the training or test set. The highest Tanimoto similarities of these compounds to their matched cholinergic or inactive (control) scaffolds were used as covariates within a random forest prediction model, with the inclusion of the molecular weight of the compound and a fragment complexity score (to account for the size and connectivity of the chemical structure). The fragment complexity was calculated using the CDK molecular properties node within KNIME. The random forest machine learning algorithm constructed 100 random classification sub-models using an information gain split criterion. The accuracy of the resulting model was tested against the test set of compounds (10% of data including actives and inactive controls) that were withheld from the initial scaffolding step. Fingerprint-based prediction models were developed as described above for the scaffold-based model with inclusion of the structural MACCS fingerprint (described above) of the molecule (for the fingerprint-based model) as the sole covariate within the random forest algorithm.The final hybrid scaffold/fingerprint-based model was constructed using composite predictions from the scaffold- and fingerprint-based models as sole covariates. Hybrid models incorporating docking binding energies for the nAChR and mAChR were constructed similarly (see below), with the inclusion of an initial clustering step where compounds were hierarchically clustered based on their binding energies. This was undertaken via calculation of Manhattan distance matrices for the binding energies of the active and inactive control compounds within the training set, with clusters assigned using a threshold distance cut-off of 0.09. Scaffolds were then calculated within 7

each binding-energy cluster and the scaffold-based random forest model constructed using the algorithms described above. The output of this model, together with that from the fingerprint-based model (detailed above) was used as the sole covariates within the final hybrid docking/scaffold/fingerprint-based model. All models were tested against the reserved (test) dataset representing 10% of the compiled compounds (both active and inactive controls) that were not used in developing the initial model. The concordance with the cholinergic status of each compound was assessed as described below in the ‘statistical analysis’ section.

Receptor-ligand docking Molecular docking of the cholinergic compounds and associated controls to the crystal structures of the nAChR and mAChR as well as AChE was undertaken using the AutoDock Vina molecular docking system (Trott and Olson, 2010) implemented in YASARA structure (ver. 17.8.15; YASARA Biosciences). The database of cholinergic compounds and associated controls were encoded as SMILES patterns with the 3D-coordinates generated and optimized using a custom YANACONDA script (available in supplemental materials) that imported each structure into YASARA and optimized its structure in vacuo via energy minimization using the NOVA force field and a simulated annealing algorithm. Following this initial round of optimization, the bond geometries were further optimized using the MOPAC implementation within YASARA. The resulting structures were stored as an SDF (structure-data file) library for docking to the respective ligand. For docking the ligands to the cholinergic targets, the structures of the human α4β2 nAChR (PDB ID 5kxi (Morales-Perez et al., 2016), human M4 mAChR (PDB ID 5dsg (Thal et al., 2016) and murine AChE in complex with Sarin (PDB ID 2whq; Sussman et al., 2009) were downloaded from the protein data bank (PDB). Following acquisition of the protein structures, each structure was protonated at pH 7.4 and visually inspected for usability as evidenced by the lack of main chain breaks and presence of well-defined residues within the ligand binding sites. The resulting structures were optimized for docking via deletion of the crystallographic waters and heteroatoms (with the exception of any bound ligands) and subjected to 8

a round of energy minimization (with the backbone atoms fixed) using the Amber14 force field along with TIP3 water molecules within a simulation cell extending 5 Å from the edge of each protein (with the addition of counter-ions (Na+ and Cl-) to neutralize the mix). Following refinement via energy minimization, simulation cells for undertaking docking analyses were setup for the assessed binding sites within each receptor, with the cell boundaries and centers detailed in Supplemental Table S3. High throughput docking of the ligand libraries to each binding site was undertaken using the ligand screening algorithm within YASARA structure, with each ligand docked to the receptor 8 times using the AutoDock Vina implementation within YASARA. The resulting poses were evaluated and the highest scoring pose reported in the final analysis. Since the nAChR has two potential ligand binding sites (the endogenous substrate (ACh) binding site and the channel pore), ligands were docked to both sites and the ligand conformation with the higher binding energy was used in the final analysis, with stronger interactions represented by numerically larger values (binding energy being defined as the additive inverse of the free energy of binding (i.e. -ΔGB)). The scoring function for binding energy in YASARA contrasts with other scoring functions that report a free energy of binding (ΔGB) where more negative results are indicative of stronger interactions.

Statistical analysis Concordance analyses between experimental cholinergic assays and predictions from our machine learning models were undertaken to assess the sensitivity, specificity and balanced accuracy of the models as well as their positive and negative predictive values. Balanced accuracy was used as a measure of the overall accuracy of the models since the active and control datasets contained differing number of substances. All data analyses were undertaken in KNIME and GraphPad Prism (ver. 6.07), with statistical significance determined at an alpha level of 0.05. For column plots, data show the mean ± SEM, with statistical significance assessed via two-tailed t-tests.

Results 9

Developing a database of cholinergic compounds and controls To capture the structural diversity of compounds that interact with the cholinergic system, on-line public databases and the published literature were curated to compile databases of cholinergic substances together with associated inactive controls from ToxCast. In a few cases, positive and control data from the ToxCast in vitro assays were not in agreement with data from other sources (e.g. the phosphothionates were not flagged as AChE inhibitors due to the in vitro assay lacking a key metabolism step needed to convert these compounds to their oxon form). These discrepancies were addressed via manual assessment of the data by subject matter experts to minimize confounders. By excluding compounds from the negative (control) dataset if they were already present within the positive dataset, the models are intentionally built to favor higher sensitivity to minimize the likelihood of “false negatives”. Thus, the models are unlikely to miss candidate molecules that have the potential to alter cholinergic neurotransmission, which can be assessed further in follow-up in vitro/in vivo assays.

Using machine learning to align cholinergic compounds with their target receptors Following scaffolding of the compiled cholinergic datasets to cover substances that interact with the three main proteins mediating cholinergic neurotransmission (nAChR, mAChR and AChE), a finalized set of 500-900 scaffolds were identified that cover the structural diversity of the cholinergic landscape. These scaffolds, together with the calculated MACCS fingerprints, were used as covariates in hybrid supervised machine-learning models to develop artificial intelligence (AI)-based classification systems to identify the likelihood an uncharacterized compound would interact with the mAChR, nAChR or AChE. When constructing the models, 10% of the data (actives and inactive controls) were randomly withheld as a test set, with the remaining 90% used to train the models. This was repeated 5 times, re-training the model with different subsets of the data withheld within the test set following which, the average concordance statistics were calculated. The scaffold/fingerprint-based models specific for the mAChR, nAChR and AChE were extremely sensitive (98.4-100%), very specific (93.7-100%) and highly accurate (96.1-100%; balanced accuracy) 10

when challenged with the reserved test set (Table 1). These models can also be used to predict whether an unknown will or will not interact with these proteins, showing both high positive (98.6-100%) and negative (93.0-100%) predictivity (Table 1). As seen in Table 1, the model for AChE was able to correctly predict the likelihood a compound would or would not interact with the enzyme with 100% concordance, highlighting the ability of the calculated scaffolds and fingerprints to identify the unique structural attributes of compounds known to inhibit AChE. The structural diversity of compounds that interact with the mAChR or nAChR was similarly captured by the scaffolding and fingerprinting algorithms, yielding >98% sensitivity for both receptors. These findings demonstrate the feasibility and utility of fingerprinting and structural scaffolding within a machine learning context to align previously uncharacterized/novel compounds with the major components of the cholinergic system.

Utility of protein-ligand docking in identifying cholinergic substances The ability of a substance to transiently interact with a receptor and drive downstream cellular events is driven primarily via a combination of non-covalent interactions (hydrogen bonds, salt bridges and Van der Waals interactions) that stabilize unique conformations (termed “poses”) of the ligand within its cognate binding site(s) (reviewed in (Lommerse and Taylor, 1997; Wijeyesakere and Richardson, 2017). The presented models using fingerprinting and scaffolding do not consider either the nature of the binding site(s) of the cholinergic target(s) or the mode of interaction between ligands and their target receptors, which requires three-dimensional modelling techniques such as ligand-receptor docking. Docking adds information on the mode of interaction of a test compound, allowing one to discern/visualize atomistic interactions (or lack thereof) within the receptor binding site, providing insights into the nature of the noncovalent interactions that mediate ligand binding. We assessed the utility of ligand-receptor docking in predicting the cholinergic potential of a compound. We docked the compounds in our nAChR, mAChR and AChE databases along with associated ToxCast controls to the binding sites on the respective receptors, with the resulting ligand poses assessed 11

visually. Compounds that interact with the nAChR and mAChR displayed significantly higher affinities (assessed via their binding energies) relative to control compounds for the respective target receptor (p < 0.001) (Fig. 3a and b). This is consistent with substances that specifically interact with a receptor displaying higher in vitro binding potencies relative to control compounds. Inclusion of a ligand binding energy factor within the initial scaffolding stage allowed for the construction of hybrid docking/scaffold/fingerprint-based random forest machine learning models. When challenged with the reserved test set of compounds, these models predicted the likelihood of a compound interacting with the mAChR or nAChR with virtually 100% sensitivity and ≥ 91.0% balanced accuracy (Table 2). As seen in Table 2, inclusion of the docking binding energies greatly improved the predictive power of our mAChR model (initial results in Table 1), yielding perfect (100%) concordance. In the case of the nAChR, inclusion of the docking results improved the model sensitivity and negative predictive value (99.6 and 97.8%, respectively) (Table 2), minimizing the likelihood of “false negatives” and demonstrating the confidence with which non-binders can be identified. The increased diversity of compounds that target the nAChR, coupled with the presence of multiple potential sites for ligand interaction within the receptor (Kouvatsos et al., 2016), results in the hybrid docking/scaffold/fingerprint-based model having slightly lower specificity relative to the nAChR model that does not consider ligand binding energies (Table 2 vs Table 1). However, it is important to note that docking results have the added benefit of yielding output (receptor-ligand complexes) that can be visually inspected to yield mechanistic and structural insights into the modes of action of compounds and possible adverse outcome pathways that may be propagated by the ligand-receptor interaction. In contrast to findings with the nAChR and mAChR compounds (Fig. 3a and b), AChE inhibitors did not display higher affinities towards AChE relative to the negative controls (Fig. 3c). It is noteworthy that AChE inhibitors display significantly (p< 0.001) lower affinities towards the enzyme when compared to the inactive control compounds (Fig. 3c), with the mean binding energy for the AChE inhibitors being 9.9 kcal/mol, indicative of unfavorable binding energetics (binding energy < 0 kcal/mol). This observation is likely due to the fact that a major subset of the known inhibitors [including organophosphorus esters 12

(Makhaeva et al., 2013) and carbamates (Wilson, 2010)] have a different interaction with AChE, acting via the formation of covalent adducts that derivatize the catalytic serine nucleophile. The docking algorithm used in this investigation does not model the transition state of the inhibitorprotein complex (i.e. it does not simulate the formation of a covalent bond between the electropositive center of the inhibitor and the active site nucleophile). This suggests that the use of binding energetics in identifying potential AChE inhibitors requires further investigation, likely involving the use of alternative techniques (e.g. the inclusion of quantum mechanical parameters to model the redox state of the reactants and formation of the serine-inhibitor bond). Moreover, addition of molecular dynamics to model conformational changes within the active site chamber of AChE following release of the leaving group (which could alter conformations of amino acids that form the active site chamber of the enzyme) may further improve the utility of ligand docking when investigating covalent inhibition of enzymes such as AChE. Since these techniques are beyond the scope of this investigation, and given that the lower binding energies of AChE inhibitors relative to controls (Fig. 3c) may confound the interpretation of the model results that incorporate binding energies for inhibitor-AChE complexes, we did not construct a hybrid docking/scaffold/fingerprint-based model for AChE, but retained the scaffold/fingerprint model.

Development of a pan-cholinergic model In addition to modeling the individual components of the cholinergic system (Table 1), we sought to understand the utility of developing a pan-cholinergic model to identify compounds that might interact with any of the three proteins that mediate cholinergic signaling. Using our master cholinergic and inactive control datasets, we constructed a hybrid scaffold/fingerprint-based model for the overall cholinergic system. As seen in Table 1, the pan-cholinergic model was able to identify cholinergic substances with unprecedented (98.7%) sensitivity and high accuracy (91.1%). Since this model includes diverse compound classes (nicotinics, muscarinics and AChE inhibitors), it displays a lower specificity relative to the other models. However, the high positive and negative predictive values (98.6 and 85.2%, respectively) of the 13

model demonstrates the high precision with which cholinergic compounds can be identified. While the pancholinergic model displays high sensitivity, its lower specificity relative to models for the individual protein targets demonstrates that while a multi-protein system can be modelled in its entirety, the structural diversity of compounds that interact with the system impairs model accuracy.

Impact of number of active compounds in training set on the utility of machine learning algorithms We also sought to understand how the size of the training set affects the predictive power of our models since we attempted to curate most substances known to target the cholinergic system. The master pan-cholinergic dataset was partitioned such that 1-90% of the cholinergic agents were randomly withheld (in 1% increments) within the training set together with 90% of the inactive controls. These compounds were used to build the fingerprint-based pan-cholinergic model, with the remaining compounds used to challenge the model. Given the size of the cholinergic datasets and the time taken to build the hybrid scaffold/fingerprint-based models (>1 day/run), only the fingerprint-based models were evaluated in this analysis. Fingerprint models were judged to be sufficient to gain some understanding on the stability of these models as a function of dataset size. As seen in Fig. 4, only 7, 4 or 12% of actives were used to train the nAChR, mAChR and AChE models respectively. The relatively small proportion of actives needed to obtain 90% sensitivity for the mAChR model supports the limited structural diversity of substances that interact with this receptor, in good agreement with the hybrid models for this target being the most accurate (Tables 1 and 2). It is interesting that a relatively greater proportion of curated “actives” (12%) were needed for the fingerprintbased AChE model (which, when combined with identified scaffolds, yielded 100% concordance) to reach 90% sensitivity. This is likely due to this analysis only involving the fingerprint-based portion of the model. Since AChE inhibitors fall within a relatively restrictive set of chemical classes, it would be expected the addition of scaffolding information would be a better predictor than a generic fingerprint that covers a broad range or organic functional groups (which would add more “noise” to the outcome).

14

Similar to findings were observed with the pan-cholinergic model, with 90% sensitivity reached when only 821 active compounds (4% of cholinergic dataset) were used to train the model (Fig. 4a).

Impact of training set diversity on the utility of machine learning algorithms The findings thus far show that highly sensitive models for the entire cholinergic system or its individual components can be built with ≤12% of the compiled database of active compounds that covers the structural diversity of the cholinergic landscape (Fig. 4). To better understand the specific impact of training the generic machine learning workflow with limited data, we re-built the hybrid scaffold/fingerprint-based models using only the positive and negative in vitro data compiled from ToxCast (17, 101 or 68 actives covering 2043 substances tested for binding the nAChR, mAChR or AChE, respectively). As was the case with the other models, 10% of the dataset was withheld as a challenge set, with the remaining 90% used to train the models. The resulting concordance measurements are shown in Supplemental Table S4. As seen in Supplemental Table S4, in contrast to the results with using all available data, no scaffold/fingerprint-based model could be built for the nAChR using the ToxCast data alone as the “actives” dataset was too small. Models could be built for mAChR and AChE, albeit with reduced sensitivities relative to those built with all available data. This is likely due to the limited structural heterogeneity of active compounds within the mAChR and AChE training sets. These findings demonstrate the limitations of using machine-learning models in cases where relatively low numbers of structurally heterogeneous active compounds are available/compiled for a target of interest. This suggests that a dataset requires sufficient representation of active and inactive (control) structural scaffolds for these machine-learning modelling techniques to reliably detect a wide array of active substances.

Use of the cholinergic profilers to screen a microelectrode array (MEA) database Recently, the microelectrode array (MEA) assay has been proposed as an in vitro technique to study changes in neurosignalling following exposure to xenobiotics. To assess the utility of implementing our 15

cholinergic profilers, we screened a recently published database of >1000 compounds from ToxCast phase II with MEA data (detailed in Strickland et al. (2018) using the scaffold/fingerprint models for the nAChR, mAChR and AChE, with the results shown in Table 3. The models flagged all the nicotinic or muscarinic compounds identified by the authors of the MEA manuscript (Table 3 right-most column). Our model for AChE inhibitors identified 55 of the 57 AChE inhibitors identified by the authors. The 2 remaining compounds [Daminozide (CAS 1596-84-5) and Aldicarb oxime (CAS 1646-75-9)] were in fact likely not AChE inhibitors, but rather Daminozide is a histone deacetylase inhibitor (Rose et al., 2012), whereas Aldicarb oxime is an inactive hydrolysis product of the carbamate Aldibcarb. In addition to identifying the previously characterized cholinergic substances, our profilers flagged many additional compounds that were active in the MEA assay that could potentially interact with these cholinergic targets. Whereas only one compound [Zamifenacin (CAS 127308-82-1)] was previously flagged as mAChR active, our profiler identified an additional 121 compounds that may interact with this target. Similarly, an additional 33 and 63 substances that interact with the nAChR and AChE were flagged by the profilers. Thus, large datasets can be rapidly screened to identify previously uncharacterized potential molecular initiating events or promising lead molecules for receptors/biological targets of interest. As a second example application, by virtue of being implemented as extensible KNIME workflows, the cholinergic profilers were incorporated as statistical parameters in a mechanistic quantitative structure activity relationship model for acute oral toxicity (manuscript in preparation).

Discussion We demonstrate the feasibility of using automated workflows to compile large datasets of compounds known to interact with the cholinergic system and the value of integrating fingerprinting, scaffolding and docking in the development of highly accurate predictive models that can rapidly identify potential cholinergic substances in large datasets.

16

Given the role of the cholinergic system as a major target for therapeutics and unintentional xenobiotic exposures, these models are expected to have broad applicability across diverse fields including drug discovery, emergency medicine and product safety evaluations. When implemented as part of an integrated testing strategy, compounds flagged by these profilers can be prioritized for further downstream in vitro or in vivo testing. This value is evident when evaluating potential lead candidate molecules or novel research samples, where the ability to identify and exclude substances with unintended cholinergic effects is important to minimize an off-target of toxicity. It is important to emphasize the value of having biological targets tested in large structurally diverse in vitro testing programs. Given that a minimal proportion (≤12%) of the compiled actives were needed to yield >90% sensitivity, it is likely that the predictive power of these models will not change for the foreseeable future. To the best of our knowledge, this work represents the first study to compile large, structurally diverse datasets encompassing virtually all compounds known to target the cholinergic system. In contrast to previous modeling practices, the presented approach used data on both active and inactive control compounds, allowing for assessment of whether a test substance interacts with a receptor or not. A significant feature of these models is incorporation of negative (inactive control) data, which is particularly important when evaluating substances that are not profiled as cholinergic since the ability to confidently assign a ‘no’ is dependent on the presence of inactive control substances within the training set. By developing and implementing scaffold/fingerprint or scaffold/fingerprint/docking models as part of broadly generalizable workflows, they can be applied to any neuronal target for which enough active and inactive control data are available. Nevertheless, the need for oversight by subject matter experts trained in neurobiology is important when implementing such models, as is curating as many active compounds as possible because omission of key chemical classes would cause failure to identify certain scaffolds. An understanding of the in vitro assay methodologies used to generate the data is essential in minimizing the incorporation of “false negative” or “false positive” results (e.g. confounding results from compounds requiring metabolism or surfactants) into the training set. In instances when both positive and negative in vitro results are present for a compound and the positive results appear reasonable, we propose that 17

preference be given to the positive result to enhance model sensitivity and minimize the likelihood that potential active compounds would be missed. We propose that when modeling any target, the scaffold/fingerprint model be built along with the scaffold/fingerprint/docking model (in cases where docking is feasible). The scaffold/fingerprint models display high accuracy and can be readily implemented without access to the software or expertise needed to undertake receptor-ligand docking. Inclusion of docking generally improves model sensitivity and allows for visualization of the ligand-receptor complexes, which can be valuable in the screening of lead molecules for the receptors of interest. Whereas this investigation focused on cholinergic compounds, the described modeling approaches are broadly generalizable and expected to provide similar benefits (predictive framework with mode of action/adverse outcome information) for any biological target with in vitro data to allow for sufficient coverage of the diversity of the chemical landscape of actives and inactive controls. Indeed, in developing models for other (unrelated targets), we have found that the accuracy of machine learning models can be limited by a range of factors such as: 1) limited structural diversity of actives (as shown in this investigation) or controls (e.g. in cases where the receptor has not been tested in an in vitro high throughput screening assay); and 2) promiscuity of the target(s) being modelled (e.g. the cholinergic inhibitors comprises a defined subset of chemical space in contrast to other targets in contrast to substances that can inhibit more complex organelles such as the mitochondria (manuscript in press). Thus, the amount of data needed to train a usable model will depend on the target(s) being modelled and should be evaluated to ensure sufficient coverage of the chemical space of known binders and non-binders (controls). Towards this end, these techniques have been readily adapted to several other important biological targets with similar results (manuscripts in preparation).

Conflict of interest The authors declare no conflicts of interest.

18

References

Bento, A.P., Gaulton, A., Hersey, A., Bellis, L.J., Chambers, J., Davies, M., Krüger, F.A., Light, Y., Mak, L., McGlinchey, S., Nowotka, M., Papadatos, G., Santos, R., Overington, J.P., 2014. The ChEMBL bioactivity database: an update. Nucleic Acids Research 42, D1083D1090. Berthold, M.R., Cebron, N., Dill, F., Gabriel, T.R., Kötter, T., Meinl, T., Ohl, P., Sieb, C., Thiel, K., Wiswedel, B., 2008. KNIME: The Konstanz Information Miner. 319-326. Bikadi, Z., Simonyi, M., 2003. Muscarinic and nicotinic cholinergic agonists: structural analogies and discrepancies. Curr Med Chem 10, 2611-2620. Bird, M.J., Thorburn, D.R., Frazier, A.E., 2014. Modelling biochemical features of mitochondrial neuropathology. Biochim Biophys Acta 1840, 1380-1392. Chen, X., Liu, M., Gilson, M.K., 2001. BindingDB: a web-accessible molecular recognition database. Comb Chem High Throughput Screen 4, 719-725. Dix, D.J., Houck, K.A., Martin, M.T., Richard, A.M., Setzer, R.W., Kavlock, R.J., 2007. The ToxCast program for prioritizing toxicity testing of environmental chemicals. Toxicol Sci 95, 5-12. Durant, J.L., Leland, B.A., Henry, D.R., Nourse, J.G., 2002. Reoptimization of MDL keys for use in drug discovery. J Chem Inf Comp Sci 42, 1273-1280. Fasoli, F., Gotti, C., 2015. Structure of neuronal nicotinic receptors. Curr Top Behav Neurosci 23, 1-17. Gaulton, A., Bellis, L.J., Bento, A.P., Chambers, J., Davies, M., Hersey, A., Light, Y., McGlinchey, S., Michalovich, D., Al-Lazikani, B., Overington, J.P., 2012. ChEMBL: a large-scale bioactivity database for drug discovery. Nucleic Acids Res 40, D1100-1107. Gilson, M.K., Liu, T., Baitaluk, M., Nicola, G., Hwang, L., Chong, J., 2016. BindingDB in 2015: A public database for medicinal chemistry, computational chemistry and systems pharmacology. Nucleic Acids Res 44, D1045-1053. Harbauer, A.B., 2017. Mitochondrial health maintenance in axons. Biochem Soc Trans 45, 10451052. Kouvatsos, N., Giastas, P., Chroni-Tzartou, D., Poulopoulou, C., Tzartos, S.J., 2016. Crystal structure of a human neuronal nAChR extracellular domain in pentameric assembly: Ligand-bound alpha2 homopentamer. Proc Natl Acad Sci U S A 113, 9635-9640. Kruse, A.C., Kobilka, B.K., Gautam, D., Sexton, P.M., Christopoulos, A., Wess, J., 2014. Muscarinic acetylcholine receptors: novel opportunities for drug development. Nat Rev Drug Discov 13, 549-560. Lee, S., Barron, M.G., 2018. 3D-QSAR study of steroidal and azaheterocyclic human aromatase inhibitors using quantitative profile of protein–ligand interactions. Journal of Cheminformatics 10. Legay, C., Mei, L., 2017. Moving forward with the neuromuscular junction. J Neurochem 142 Suppl 2, 59-63. Lommerse, J.P., Taylor, R., 1997. Characterising non-covalent interactions with the Cambridge Structural Database. J Enzyme Inhib 11, 223-243. Makhaeva, G.F., Radchenko, E.V., Palyulin, V.A., Rudakova, E.V., Aksinenko, A.Y., Sokolov, V.B., Zefirov, N.S., Richardson, R.J., 2013. Organophosphorus compound esterase profiles as predictors of therapeutic and toxic effects. Chem Biol Interact 203, 231-237. 19

Morales-Perez, C.L., Noviello, C.M., Hibbs, R.E., 2016. X-ray structure of the human α4β2 nicotinic receptor. Nature 538, 411-415. Nelms, M.D., Mellor, C.L., Cronin, M.T., Madden, J.C., Enoch, S.J., 2015. Development of an in Silico Profiler for Mitochondrial Toxicity. Chem Res Toxicol 28, 1891-1902. Richardson, R.J., Hein, N.D., Wijeyesakere, S.J., Fink, J.K., Makhaeva, G.F., 2013. Neuropathy target esterase (NTE): overview and future. Chem Biol Interact 203, 238-244. Rose, N.R., Woon, E.C.Y., Tumber, A., Walport, L.J., Chowdhury, R., Li, X.S., King, O.N.F., Lejeune, C., Ng, S.S., Krojer, T., Chan, M.C., Rydzik, A.M., Hopkinson, R.J., Che, K.H., Daniel, M., Strain-Damerell, C., Gileadi, C., Kochan, G., Leung, I.K.H., Dunford, J., Yeoh, K.K., Ratcliffe, P.J., Burgess-Brown, N., von Delft, F., Muller, S., Marsden, B., Brennan, P.E., McDonough, M.A., Oppermann, U., Klose, R.J., Schofield, C.J., Kawamura, A., 2012. Plant Growth Regulator Daminozide Is a Selective Inhibitor of Human KDM2/7 Histone Demethylases. J Med Chem 55, 6639-6643. Silman, I., Sussman, J.L., 2017. Recent developments in structural studies on acetylcholinesterase. J Neurochem 142 Suppl 2, 19-25. Sterling, T., Irwin, J.J., 2015. ZINC 15--Ligand Discovery for Everyone. J Chem Inf Model 55, 2324-2337. Strickland, J.D., Martin, M.T., Richard, A.M., Houck, K.A., Shafer, T.J., 2018. Screening the ToxCast phase II libraries for alterations in network function using cortical neurons grown on multi-well microelectrode array (mwMEA) plates. Arch Toxicol 92, 487-500. Sussman, J.L., Ekström, F., Hörnberg, A., Artursson, E., Hammarström, L.-G., Schneider, G., Pang, Y.-P., 2009. Structure of HI-6•Sarin-Acetylcholinesterase Determined by X-Ray Crystallography and Molecular Dynamics Simulation: Reactivator Mechanism and Design. PLoS ONE 4, e5957. Szinicz, L., 2005. History of chemical and biological warfare agents. Toxicology 214, 167-181. Tan, Y.-J., Ting, A.E., 2000. Non-ionic detergent affects the conformation of a functionally active mutant of Bcl-XL. Protein Engineering, Design and Selection 13, 887-892. Thal, D.M., Sun, B., Feng, D., Nawaratne, V., Leach, K., Felder, C.C., Bures, M.G., Evans, D.A., Weis, W.I., Bachhawat, P., Kobilka, T.S., Sexton, P.M., Kobilka, B.K., Christopoulos, A., 2016. Crystal structures of the M1 and M4 muscarinic acetylcholine receptors. Nature 531, 335-340. Tiwari, P., Dwivedi, S., Singh, M.P., Mishra, R., Chandy, A., 2013. Basic and modern concepts on cholinergic receptor: A review. Asian Pacific Journal of Tropical Disease 3, 413-420. Trott, O., Olson, A.J., 2010. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem 31, 455-461. Verma, S., Kumar, A., Tripathi, T., Kumar, A., 2018. Muscarinic and nicotinic acetylcholine receptor agonists: current scenario in Alzheimer's disease therapy. J Pharm Pharmacol 70, 985-993. Wijeyesakere, S.J., Richardson, R.J., 2017. CHAPTER 7. In silico Chemical–Protein Docking and Molecular Dynamics. 174-190. Wilson, B.W., 2010. Cholinesterases. 1457-1478.

20

Acknowledgments: The authors thank Amanda Parks and Tyler Auernhammer for their assistance in curating datasets and helpful discussions. We also thank Drs. Johnson Thomas and Amanda Andrus for their technical review of the manuscript and helpful conversations as well as Dr. Elmar Krieger for his assistance with the development of custom YANACONDA macros.

Author contributions: SJW and DW contributed to conception of the hypothesis, design of the experiments and analysis of the data. SJW setup and ran the computational analyses to collect data to build the predictive models. SJW, DW and SM contributed to the writing of the manuscript.

21

Figure Legends Figure 1. Structures of model cholinergic ligands. Figure shows the structure of acetylcholine (ACh), the endogenous ligand/neurotransmitter that mediates cholinergic neurotransmission. The structures of nicotine and muscarine (that are known to bind the nicotinic (nAChR) and muscarinic (mAChR) ACh receptors respectively), together with sarin (a model AChE inhibitor) are shown. Figure 2. Overview of the custom KNIME algorithm used to scaffold compounds. Figure 3. Target receptors with known structures allow for the use of protein docking in model building. When interactions between a ligand and a receptor are primarily driven by non-covalent interactions, compounds that are known to interact with the receptor display a higher in silico affinity towards the receptor. Figure shows the results of docking the compiled (a) mAChR, (b) nAChR or (c) AChE binders and associated control compounds to the ligand binding sites on the indicated protein. Following docking, the binding energies for the ligand-receptor interaction were calculated and tabulated, with each panel depicting the mean binding energy for the ligand-receptor interaction ± SEM. Statistical significance was assessed via two-tailed t-tests. The affinity of the interaction between the ligands and target receptor is defined by the binding energy with larger values representing a more favorable interaction. Docking also allows visualization into the mode of interaction of the ligands with target receptors. Figure 4. Performance of the supervised machine-learning models is dependent upon the size of the training set. Plots depict the change in the sensitivity as a function of the fraction of actives within the master pan-cholinergic (A), nAChR binding (B), mAChR binding (C) and AChE inhibitor (D) datasets used in training the fingerprint-based random forest machine learning models. Models were trained with the indicated proportion of active compounds together with 90% of the controls and challenged with the remaining active and control compounds. Plots depict the mean concordance measurement (%) ± SEM for three independent assessments. Fingerprint-based models (without scaffolding) were judged to be sufficient to glean the relevance of dataset size/diversity on model performance. 22

Table 1. Concordance analysis for prediction of the cholinergic potential of compounds via the use of hybrid fingerprint/scaffold-based machine learning models that predict the ability of compounds to interact with either of the three individual proteins that comprise the cholinergic system [the muscarinic (mAChR) or nicotinic (nAChR) acetylcholine receptors or acetylcholinesterase (AChE)] or with any of the cholinergic proteins (pan-cholinergic model). Data represent the average of 5 independent runs where a random set of compounds comprising 10% of the actives and controls within each dataset was withheld for assessing each model’s performance.

Concordance (%)

mAChR

nAChR

AChE

Pan-cholinergic

Sensitivity

99.3

98.4

100

98.7

Specificity

94.3

93.7

100

83.5

Predictive Positive Value

99.4

98.6

100

98.6

Predictive Negative Value

93.4

93.0

100

85.2

Balanced Accuracy

96.8

96.1

100

91.1

Number of actives

4,481

8,931

7,400

20,536

23

Table 2. Concordance analysis for prediction of the likelihood a compound will interact with the acetylcholine receptors (nAChR or mAChR) using a hybrid docking/scaffold/fingerprint-based machine learning model. This approach involves consideration of ligand binding energies within the scaffolding algorithm. Acetylcholinesterase (AChE) was excluded from this analysis since a major subset of compounds that inhibit this enzyme require the formation of a covalent adduct and therefore, could not be docked using the available model. Data represent the average of 5 independent runs where a random set of compounds comprising 10% of the actives and controls within each dataset was withheld for assessing each model.

Concordance (%)

mAChR

nAChR

Sensitivity

100

99.6

Specificity

100

82.4

Predictive Positive Value

100

96.1

Predictive Negative Value

100

97.8

Balanced Accuracy

100

91.0

24

Table 3. Application of cholinergic profilers to a large dataset of 1,024 compounds with microelectrode array (MEA) data (Strickland et al., 2018). Table shows the number of compounds assigned by the authors of the MEA analysis as well as the number of compounds identified by the profilers described in this manuscript.

Number of compounds Assigned to cholinergic mechanism in Strickland et al.

No. of cholinergic compounds identified by Strickland et al. with positive MEA/total cholinergics in DB

Identified by cholinergic profilers

No. identified by cholinergic profilers / No. assigned cholinergic MOA by Strickland et al.

nAChR

9

2/9

33

9/9

mAChR

1

1/1

121

1/1

AChE

57

15/57

63

55/57

Target

25

Figure 1

26

Figure 2

27

Figure 3

28

Figure 4

29