Discovery of novel dual adenosine A1/A2A receptor antagonists using deep learning, pharmacophore modeling and molecular docking

Adenosine receptors (ARs) have been demonstrated to be potential therapeutic targets against Parkinson’s disease (PD). In the present study, we describe a multistage virtual screening approach that identifies dual adenosine A1 and A2A receptor antagonists using deep learning, pharmacophore models, and molecular docking methods. Nineteen hits from the ChemDiv library containing 1,178,506 compounds were selected and further tested by in vitro assays (cAMP functional assay and radioligand binding assay); of these hits, two compounds (C8 and C9) with 1,2,4-triazole scaffolds possessing the most potent binding affinity and antagonistic activity for A1/A2A ARs at the nanomolar level (pKi of 7.16–7.49 and pIC50 of 6.31–6.78) were identified. Further molecular dynamics (MD) simulations suggested similarly strong binding interactions of the complexes between the A1/A2A ARs and two compounds (C8 and C9). Notably, the 1,2,4-triazole derivatives (compounds C8 and C9) were identified as the most potent dual A1/A2A AR antagonists in our study and could serve as a basis for further development. The effective multistage screening approach developed in this study can be utilized to identify potent ligands for other drug targets.


Introduction
Parkinson's disease (PD) is a common and complex neurodegenerative disorder that is characterized by the early prominent death of dopaminergic neurons in the substantia nigra pars compacta and the abnormal aggregation of the α-synuclein protein, called Lewy bodies and Lewy neurites [1]. Decreased dopamine in the basal ganglia can cause classical motor symptoms [2], including bradykinesia, resting tremors, and postural instability [3], and nonmotor symptoms, including constipation, depression, sleep disturbance, apathy, hallucinations and dementia [4]. L-3,4-Dihydroxyphenylalanine (L-dopa), which is the direct precursor of dopamine, is commonly used as dopamine replacement therapy in the treatment of PD motor symptoms [5]. Although L-dopa can help relieve PD motor symptoms, its chronic use may cause side effects such as motor complications (motor fluctuations and dyskinesias) [6]. In addition, dopamine receptor agonists, catechol O-methyltransferase inhibitors, monoamine oxidase B (MAOB) inhibitors, amantadine and anticholinergic drugs are available on the market for the treatment of PD [7,8]. However, these drugs are mainly used to replace the concentration and/or effect of dopamine in the brain and only solve the motor symptoms but not the nonmotor symptoms [9]. Therefore, it is necessary to develop effective therapies that simultaneously solve the motor symptoms and nonmotor symptoms in medical treatments of PD.
In recent years, adenosine receptor (AR) antagonists have attracted much attention in the development of nondopaminergic therapies for the treatment of PD. Interestingly, high-resolution structures of A 1 and A 2A adenosine receptors are available [10][11][12][13][14][15][16][17][18][19][20], providing an excellent opportunity for structure-based drug design (SBDD). Thus, virtual screening efforts have identified potent antagonists [21][22][23]. Among the four human adenosine receptors (A 1 , A 2A , A 2B and A 3 ) [24], the A 2A receptor, which is a potential target for the treatment of PD, has been intensively studied. A 2A AR antagonists have been developed as a potential class of nondopaminergic antiparkinsonian agents that can relieve patients of symptoms of depression (nonmotor symptom of PD) [25]. The A 2A AR antagonist KW-6002 has been proven to have antidepressant activity in the forced swim test and tail suspension test in rodents [26,27]. In addition, epidemiological and experimental data indicate that adenosine A 2A receptor antagonists have neuroprotective effects [28]. Compared with A 2A AR, A 1 AR is more widely distributed in the central nervous system. A 1 AR is highly expressed in brain regions, including the hippocampus and prefrontal cortex. These regions are important for emotion and cognitive function. Furthermore, since A 1 AR antagonists can increase the release of acetylcholine and glutamate, and improve cognitive dysfunction [29], A 1 AR antagonism may improve the cognitive deficits experienced in PD as illustrated in animals studies [30].
Among the adenosine receptor antagonists, caffeine, which is a xanthine derivative, has been found to be a nonselective A 1 AR and A 2A AR antagonist. Epidemiological research has shown that a correlation exists between the intake of coffee or caffeine and a reduced risk of PD [31]. Subsequently, several dual A 1 /A 2A AR antagonists were further developed and proven to not only improve dyskinesia but also enhance cognition and play a neuroprotective role through the antagonistic effect of A 2A AR [3,32,33]. Dual A 1 /A 2A AR antagonists may display synergistic motor activation. A 1 AR antagonism promotes presynaptic dopamine release, while A 2A AR antagonism promotes postsynaptic dopamine release [34]. In conclusion, dual A 1 /A 2A AR antagonists not only treat the motor symptoms of PD and have neuroprotective effects but may also improve nonmotor symptoms [35]. In 2007, Mihara et al. discovered the dual A 1 /A 2A AR antagonist 5-[5-amino-3-(4-fluorophenyl)pyrazin-2-yl]-1-isopropylpyridine-2(1H)-one (ASP5854) (Fig 1A), which showed high affinity for human A 1 and A 2A receptors with K i values of 9.03 nM and 1.76 nM, respectively [36]. In 2008, ASP5854 was further confirmed to have activity comparable to that of existing anti-Parkinson's disease drugs, but no further research investigating its safety was conducted [37]. In 2010, Shook et al. designed and synthesized an arylindenopyrimidine (2-amino-4-phenyl-8-(pyrrolidin-1-ylmethyl)-5H-indeno [1,2-d]pyrimidin-5-one, Fig 1B) derivative as a dual A 1 /A 2A AR antagonist that has excellent activity in various animal models of PD after oral administration, but further studies revealed that it was genotoxic [3]. In 2015, Robinson et al. characterized a small set of 2-aminopyridines as dual A 1 /A 2A AR antagonists. These compounds are considered useful for motor diseases, such as PD. The most potent compound found from this campaign was 4-(5-methylfuran-2-yl)-6-[3-(piperidine-1-carbonyl)phenyl]pyrimidin-2-amine ( Fig 1C), which was able to bind A 1 and A 2A receptors with K i values of 9.54 nM and 6.34 nM, respectively [38,39]. Therefore, the identification of novel dual A 1 /A 2A AR antagonists is of great significance for the development of novel agents for the treatment of PD [40].
Deep learning (DL), or deep neural networks, has been developed based on artificial neural networks (ANNs), which are inspired by the biological structure and function of the brain. To date, several DL frameworks, such as deep neural networks (DNNs), convolutional neural networks (CNNs), deep belief networks (DBNs), and recurrent neural networks (RNNs), have been developed. Currently, the deep learning procedure, which is the dominant component of machine learning (ML) in artificial intelligence, has emerged as a vital tool for expediting the in silico screening of potential hits [41]. Unterthiner et al. showed that the efficacy of DL in virtual drug screening was better than that of seven other screening methods using ChEMBL benchmark data [42]. Lenselink et al. demonstrated that the screening performance of DNNs outperformed that of other machine learning methods such as support vector machines, random forests, naïve Bayes and logistic regression models [43]. Bilsland et al. trained an ANN classification model to screen for senescence-inducing compounds from a library of~2 M lead-like compounds and identified a benzimidazolone compound with a low micromolar IC 50 using in vitro assays as a potential hit for the development of selective cell cycle inhibitors [44]. Wallach et al. used CNN to design a structure-based model for the prediction of the bioactivity of small molecules for drug development and successfully predicted new active molecules for targets without known modulators [45]. Recently, Rifaioglu et al. developed DEEPScreen, which is a novel DTI prediction system based on deep convolutional neural networks that can predict new interactions between the drug cladribine and JAK proteins as confirmed by in vitro experiments using cancer cells [46]. Overall, deep learning plays an increasingly important role in modern drug discovery.
To date, the discovery of dual A 1 /A 2A AR antagonists has been achieved through structural modification of specific scaffolds to obtain new compounds, while a multistage virtual screening approach combining ligand-based and structural-based methods has rarely been used. Although the rapid development of the pharmacophore [47] and molecular docking [48] methods enabled the screening of very large databases containing billions of diverse molecules, these methods are individually unable to achieve perfect effects. For example, the inaccurate scoring functions implemented in docking methods often lead to a low hit rate and a high false positive rate. In this study, a multistage virtual screening system consisting of deep learning, pharmacophore models, and molecular docking methods was used to identify novel dual A 1 / A 2A AR antagonists from the ChemDiv database. The ChemDiv collection consists of more than 14,000 chemical families (chemotypes) and 1,250,000 diverse drug-like small molecules.
Nineteen compounds were selected and tested to determine their A 1 AR and A 2A AR functional activities; of these compounds, five compounds were found to have dual A 1 /A 2A AR antagonistic activity through cAMP functional assays and radioligand binding assays. Among the five compounds found to have dual A 1 /A 2A AR antagonistic activity, compounds C8 and C9, which contained a 1,2,4-triazole scaffold, showed the highest antagonistic activity, reaching nanomolar inhibition.

Construction and validation of the DNN and CNN models
In this study, we used the self-developed Python script based on the algorithm from reference 50 to construct the DNN and CNN classification models of the dual A 1 /A 2A AR antagonists using a library containing 310 bioactive dual A 1 /A 2A antagonists (K i < 40 nM) and 405 nonbioactive antagonists (K i > 1000 nM). The models used the extended connectivity fingerprint 4 (ECFP4) [49] and neural fingerprint (NFP) [50]. Six evaluation indicators, including sensitivity (SE), specificity (SP), prediction accuracy of active molecules (Q+), prediction accuracy of inactive molecules (Q-), Matthews correlation coefficient (MCC) and area under the curve (AUC), were used to evaluate the classification ability of the models.
To assess the impact of the different batch sizes on the performance of the DNN and CNN classification models, we constructed the DNN and CNN models using batch sizes from 50 to 300 with an interval of 50. Table 1 lists the statistical evaluation results of the DNN classification models based on the ECFP4 of the test set under different batch sizes. All DNN classification models show very good SE, SP, Q+, Q-, MCC and AUC values. Among the six DNN classification models, Model_D5 based on ECFP4 has the best prediction ability, with an MCC of 0.891 and an AUC of 0.997. The optimized hyperparameters of Model_D5 which had the best performance, are as follows: batch size = 250, learning rate = 0.001, num_epochs (number of training epochs) = 500, dropout = 0.2, L2_regulation_type (L2 regularization parameters) = 0.0001, three hidden layers with layer widths of 3000, 2000 and 1000, and a final fully connected layer that uses the Softmax algorithm to produce classification results.
We further analyzed the performance of the CNN models in predicting dual A 1 /A 2A AR antagonists. The statistical evaluation of six CNN models based on NFP is shown in Table 2. The table shows that the batch size has a great impact on the prediction ability of the CNN classification models. Model_C6, which has a batch size of 300, shows the lowest MCC and AUC values, while Model _C4, which has a batch size of 200, shows the highest MCC and AUC values. The optimized hyperparameters of Model_C4, which has the best performance are as follows: learning rate = 0.001, num_epochs (number of training epochs) = 500, batch size = 200, L2_regulation_type (L2 regularization parameters) = 0.0001, and fingerprint_network_architecture (convolutional layers) = 5. Based on the evaluation results of the DNN and CNN classification models, Model_D5 and Model_C4, which had the best performance, were selected, and we conducted the following virtual screening of dual A 1 /A 2A AR antagonists.

Pharmacophore model generation and validation
A training set containing 14 active ingredients (S1 Table) was used to generate the dual A 1 / A 2A AR antagonist pharmacophore models. Eleven hypotheses were generated, which matched 9-13 of the 14 activities. To evaluate pharmacophore hypotheses, a validation set of 42 dual antagonists of the A 1 AR and A 2A AR subtypes and 913 decoy compounds was used to explore the ability of the pharmacophore hypothesis to distinguish the dual A 1 /A 2A AR antagonists from the decoys. The statistical results of predicting the dual antagonists in the validation set are summarized in Table 3. The EF1% (enrichment factor in the top 1% of compounds screened) and BEDROC (α-160.9) (Boltzmann-enhanced discrimination of receiver operating characteristic) were used as "early recognition" metrics [51]. As shown in Fig 2, the AADR_4 and AAADR_1 hypotheses shared the same four pharmacophore sites including two hydrogen bond acceptors (A), one hydrogen donor (D) and one aromatic ring (R); however, the AAADR_1 hypothesis contains one additional hydrogen bond acceptor at the far end. When a minimum of four sites of the AAADR_1 model were used for screening (including the AADR_4 model), the most dual A 1 /A 2A AR antagonists (39 of 42) were retrieved from the validation set, corresponding to an ROC of 0.89. Finally, the five-pharmacophore-site (AAADR_1) hypotheses, which matched 12 of the 14 active ingredients in the training set, were selected for further virtual screening.

Virtual screening
Our previous study showed that a multistage approach sequentially integrating machine learning classification models and pharmacophore and molecular docking methods is efficient in identifying bioactive compounds from chemical databases [52]. In this study, the performance of the hierarchical multistage virtual screening approach was evaluated through a validation set containing 433 bioactive dual A 1 /A 2A AR antagonists (40 nM < K i <600 nM) and 11605 decoys (S2 Table). Then, DNN and CNN classification models, a pharmacophore model (AAADR_1) and molecular docking were applied for the discovery of dual A 1 /A 2A AR antagonists by performing multistage virtual screening against the ChemDiv library (1,178,506 compounds) (Fig 3). In the first stage, the DNN and CNN classification models were separately utilized to filter 1,178,506 compounds. Initially, 209,585 and 171,233 compounds passed the filter that applied the DNN and CNN classification models, respectively. In total, 58,886 compounds simultaneously predicted by the DNN and CNN classification models were retained. In the second stage, a pharmacophore model was applied to remove the compounds that were not adequate for a minimum of four sites, and 5,629 of the 58,886 compounds that matched the requirements of the AAADR_1 model were retained. In the third stage, the 5,629 compounds were subjected to molecular docking screening based on the X-ray structures of A 1 AR (PDB ID 5N2S) and A 2A AR (PDB ID 5IU4) using Glide HTVS, SP and XP functions. After the XP molecular docking, 43 compounds that simultaneously bind A 1

Biological evaluation
Functional experiments. Initially, the biological activity of selected compounds C1 -C19 towards A 1 and A 2A ARs was evaluated by a cAMP functional assay. The potency of the antagonists at the A 1 and A 2A ARs is depicted in Table 4, which shows that eight compounds (C4, C7, C8, C9, C10, C15, C17 and C19) of the 19 compounds (Fig 4) display potent antagonistic  (Fig 5B). The concentration-response curves of the eight compounds in the radioligand binding assay are depicted in S4 and S5 Figs. Table 4 shows that compounds C8 and C9 also possess the highest binding affinity for A 1 /A 2A ARs with K i values at the nanomolar level (pK i values in the range 7. 16-7.49), which is consistent with the highest antagonistic potency with pIC 50 values in the range 6.31-6.78. Overall, these experiments confirm that compounds C8 and C9 are the most potent dual A 1 /A 2A AR antagonists among our selected compounds.

Compound number
Compound ID cAMP assays (pIC 50 ) Binding affinities (pK i ) Pharmacophore matching

Novelty of the new hits
To evaluate the novelty of the identified dual A 1 /A 2A AR antagonists, the pairwise similarity between compounds C8 and C9 and known A 1 /A 2A AR antagonists, including all dual antagonists found in ChEMBL25 [53], was calculated using Morgan, ECFP4 and MACCS fingerprints in KNIME [54]. The molecular similarity between the two molecules was compared

Molecular modeling exploration
To explore the dual antagonistic activities of compound C8, C9 exhibited nanomolar inhibition against both A 1 AR and A 2A AR, and molecular docking using Glide and MD simulation studies were carried out. In our docking studies, C8 and C9 had similar binding energies against A 1 AR and A 2A AR (Table 5). In total, 4 complexes of C8 and C9 bound to A 1 AR and A 2A AR were generated, and each complex was embedded in a hydrated phospholipid bilayer including POPE lipid molecules and subjected to 100-ns MD simulation using the AMBER ff14sb force field. After 100 ns of the MD simulation we obtained a stable bilayer ( The stability of the ligand inside the binding area was assessed by measuring its RMSD (RMSD lig ). Among the four docking poses used as starting structures, we obtained RMSD lig less than 2 Å (Table 5 and S9 Fig), suggesting stable binding in agreement with the experimental binding and functional data. In addition, we calculated the binding free energy of C8 and C9 with A 1 AR and A 2A AR by using the MM-GBSA method based on the last 5 ns trajectory of the MD simulations ( Table 5). The calculated binding free energy of the C8-A 1 AR complex (ΔG eff = -41.67 kcal/ mol) is similar to that of the C8-A 2A AR complex (ΔG eff = -42.22 kcal/mol). In addition, the calculated binding energy of the C9-A 1 AR complex (ΔG eff = -43.55 kcal/mol) is similar to that of the C9-A 2A AR complex (ΔG eff = -40.70 kcal/mol). Therefore, the binding free energies show that compounds C8 and C9 have similar affinities to A 1 AR and A 2A AR, which is consistent with the results of the functional assay and radioligand binding assays.
Binding modes exploration. The 100-ns MD simulations showed that compounds C8 and C9 were stabilized with similar orientations in the binding pockets of both A 1 AR and A 2A AR (Fig 6). In the C8-A 1 AR complex and the C9-A 1 AR complex (Fig 6A and 6B), the 5-amino group and the nitrogen of the 1,2,4-triazol were hydrogen bonded to the side chain carbonyl and amino group of N254 6.55 , respectively. The 1,2,4-triazol formed π-π stacking interactions with the phenyl group of F171 ECL2 . These interactions are consistent with the observations in the mutagenesis experiments that N254 6.55 and F171 ECL2 play important roles in antagonist binding. These interactions were maintained from the origin to the end of the 100-ns simulation (S10 Fig). Additionally, the 3-methoxyphenyl (or 3-chorophenyl) exhibited nonpolar interactions with the side chain of I274 7.39 , and the 2-chorophenyl extended deeper into the other side of the pocket formed by the hydrophobic side chains of V87 3.32 , L88 3.33 and W247 6.48 . In the C8-A 2A AR complex and the C9-A 2A AR complex (Fig 6C and 6D), the 5-amino group and the oxygen of the methanone formed hydrogen bonds with the side chain carbonyl and amino group of N253 6.55 . Additionally, the 5-amino group formed a hydrogen bond with the carboxyl side chain of E169 ECL2 . The 1,2,4-triazol interacted with the phenyl group of F168 ECL2 through π-π stacking. These interactions were maintained throughout the 100-ns simulation (S11 Fig). The 2-chlorophenyl formed hydrophobic interactions with Y271 7.36 , and the 3-methoxyphenyl (or 3-chorophenyl) extended deeper into the orthosteric binding area and formed hydrophobic interactions with residues V84 3.32 , L85 3.33 , L249 6.51 and I274 7.39 . Although the orientation of C8 (or C9) was different in the binding pocket of A 1 AR and A 2A AR, the number of hydrogen bonds, π-π stacking and hydrophobic interactions were To compare the binding differences between A 1 AR and A 2A AR of nondual A 1 /A 2A AR antagonists, we also performed MD simulations against compound C10, which has stronger binding affinity and antagonistic activity to A 1 AR than A 2A AR. In the C10-A 1 AR complex, the N-propyl group and the nitrogen of the pyrazolo[1,5-a]pyrimidinformed two hydrogen bonds with the side chain carbonyl and amino group of N254 6.55 . Pyrazolo[1,5-a]pyrimidin interacted with the phenyl side chain of F171 ECL2 through π-π stacking. In the C10-A 2A AR complex, C10 possessed a different orientation and formed only one hydrogen bond with the amino group of N253 6.55 via the nitrogen of pyrazolo[1,5-a]pyrimidin. Pyrazolo[1,5-a]pyrimidin formed π-π stacking interactions with the phenyl side chain of F168 ECL2 (Figs 7 and S12). Compared to the dual antagonists (which had similar a number of hydrogen bonds, π-π stacking and hydrophobic interactions in the two receptors), C10 had one less hydrogen bond in the C10-A 2A AR complex than the C10-A 1 AR complex, which may be the reason for its weaker binding affinity and antagonistic activity for A 2A AR.

Conclusions
In the quest for novel dual A 1 /A 2A AR antagonists as putative agents for the treatment of Parkinson's disease, we applied a multistage virtual screening approach combining deep learning, pharmacophore and molecular docking methods to screen the ChemDiv library (1,178,506 compounds) and tested 19 hits by in vitro assays. Initially, the cAMP functional assay identified five compounds that possessed antagonist activity towards A 1 /A 2A AR with pIC 50 values of 4.20-6.78. The radioligand binding assays confirmed that six of the eight compounds with antagonistic activity for A 1 AR and A 2A AR (pIC 50 of 4.20-6.78) also had consistent binding affinity against A 1 /A 2A ARs (pK i of 4.71-7.49). In particular, compounds C8 and C9 showed the highest binding affinity and functional activity for A 1 /A 2A ARs with K i values at the nanomolar level (pK i of 7.16-7.49 and pIC 50 of 6.31-6.78). Compounds C8 and C9 are novel 1,2,4-triazole derivatives acting as dual A 1 /A 2A AR antagonists. The MD simulations of the

PLOS COMPUTATIONAL BIOLOGY
Discovery of novel dual adenosine A 1 /A 2A receptor antagonists complexes between the A 1 /A 2A ARs and C8 and C9 further suggest strong binding interactions. Therefore, compounds C8 and C9 are hits with potential after optimization for the development of anti-Parkinson's disease agents.

Data set preparation
Dual antagonists with binding constants (K i ) for A 1 AR and A 2A AR subtypes were investigated in the ChEMBL25 database [53]. A dataset consisting of 310 bioactive dual antagonists from approximately 100 series (K i < 40 nM) and 405 nonbioactive compounds in the same series (K i > 1000 nM) were used to train the DNN and CNN classification models. The dataset was divided into a training set, test set and validation set at a ratio of 7:2:1.
Canvas similarity and clustering from Schrödinger were used to calculate the Tanimoto coefficients (Tc) of similarity between every pair of structures among 310 dual antagonists (K i < 40 nM). Based on the calculated Tc values, similar structures were grouped together into the same cluster, yielding 14 clusters. Then, 14 representative structures of the diverse clusters were subsequently used to build dual antagonist pharmacophore models of the A 1 AR and A 2A AR subtypes. To evaluate the performance of the pharmacophore models, an extra validation set comprising 42 dual antagonists (K i < 40 nM) of A 1 AR and A 2A AR subtypes and 913 decoy compounds was applied. With respect to known dual antagonists, the decoys were selected from the ZINC database using DecoyFinder [60] based on the following criteria: (1) Tanimoto coefficient < 0.75, (2) number of hydrogen bond acceptors ± 2, (3) number of hydrogen bond donors ± 1, (4) MW ± 25, (5) logP ± 1, and (6) Tanimoto coefficient less than 0.9 with respect to the other decoys.
All compounds in the training set, test set, validation set and ChemDiv library were prepared using Schrodinger's Ligprep (version 10.2, Schrödinger, LLC) with the default settings to convert 2D structures to 3D structures, add hydrogen atoms, and generate tautomers, stereoisomers and protonation states at pH 7.0 ± 2.0 using Epik (version 4.3, Schrödinger, LLC).

DNN and CNN
A deep neural network (DNN) (Fig 8A) is a neural network with a certain complexity. DNNs comprise multiple levels of nonlinear operations and have many hidden layers. Neurons transfer information or signals to other neurons according to the received input, forming a complex network, and learning through some feedback mechanism [61]. A convolutional neural network (CNN) (Fig 8B) is a feed-forward neural network consisting of one or more convolutional layers and a fully connected layer at the top (corresponding to a classic neural network). CNNs also include associated weights and a pooling layer [62]. This structure enables the convolutional neural network to utilize the two-dimensional structure of the input data. Compared with other deep learning structures, convolutional neural networks can provide better results in image and speech recognition. Compared with other deep, feed-forward neural networks, convolutional neural networks need to consider fewer parameters, rendering these networks attractive deep learning structures.
Various statistical parameters were used to validate the performance of the DNN and CNN models. To evaluate the classification ability, parameters such as sensitivity (SE), specificity (SP), prediction accuracy of active molecules (Q+), prediction accuracy of inactive molecules (Q-), Matthews correlation coefficient (MCC) and the area under the receiver operating characteristic curve (AUC), were calculated using the following equations [63]: MCC ¼ TP � TN À FN � FP ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where TP, FP, TN and FN denote true positives, false positives, true negatives and false negatives, respectively. In this study, we aimed to improve the performance of dual A 1 /A 2A AR antagonist identification by deep learning. L2 regularization and "dropout" layers were used to prevent overfitting, and the rectified linear unit (ReLU) [64] was adopted as the activation function. The adaptive moment estimation (Adam) method [65] was used as an optimization algorithm for stochastic gradient descent in the process of training the deep learning models. In general, the performance of deep learning models is sensitive to the appropriate hyperparameters. To select the optimal hyperparameters for the deep learning models, a grid search approach [66] was used to optimize the hyperparameters based on the size of the training, including the learning rate, L2 regularization and number of hidden layers, and the size of the eigenvector in the hidden layers of the CNN model. In this stage, 209,585 and 171,233 compounds from the Chem-Div library (1,178,506 compounds) passed the filter that applied the DNN and CNN classification models, respectively. In total, 58,886 compounds that were simultaneously predicted by the DNN and CNN models were retained. The code and data sets used to train and test both the DNN and CNN models can be downloaded at https://github.com/Houshujing/ Adenosine.

Pharmacophore modeling
For the dual A 1 /A 2A AR antagonist pharmacophore generation, 14 diverse structures (K i < 40 nM) selected from 310 dual antagonists were aligned in 3D space to generate common feature pharmacophores using the Phase module (version 5.4, Schrödinger, LLC). Each antagonist structure was minimized using the OPLS-2005 force field [67] implemented in Macromodel (version 11.9, Schrödinger, LLC), and described by a list of pharmacophore sites to characterize the chemical features contributing to the compound-protein interactions between the antagonists and A 1 AR and A 2A AR subtypes. The perceived pharmacophore hypotheses were identified once the pharmacophore features were common to 50% of the antagonists. The generated pharmacophore hypotheses were examined using an extra validation set consisting of 42 dual antagonists and 913 decoy compounds. In this stage, 5,629 of the 58,886 compounds that matched the pharmacophore model were retained.

Molecular docking
A Virtual Screening Workflow (version 7.8, Schrödinger, LLC) was used to perform the molecular docking calculations using the X-ray structures of A 1 AR cocrystallized with the antagonist PSB-36 (PDB ID 5N2S) and A 2A AR cocrystallized with the antagonist ZM-241385 (PDB ID 5IU4). For each complex, the energy grid was generated at the centroid of the cocrystallized ligand. In addition, QikProp (version 5.5, Schrödinger, LLC) and Lipinski's rule were applied to eliminate molecules with undesirable drug likeness properties. The 5,629 compounds that passed QikProp and Lipinski's filters were docked and ranked sequentially using the Glide HTVS, SP and XP score functions based on the Glide score. The top 50% good scoring compounds through HTVS docking were filtered using SP. The top 25% good scoring compounds through SP docking were filtered using XP. Then, 43 of the top 25% good scoring compounds through XP against A 1 AR and A 2A AR were further screened using visual inspection.

Radioligand binding assays
The radioligand binding assays were also performed by Pharmaron (Beijing).  7.42 ). The mutant receptor residues were mutated back to the wild-type residues, the cocrystal T4 lysozyme was removed, and unnecessary small molecules were removed. The missing residues were added using the Prime module of Schrodinger. Hydrogen atoms were added through Schrodinger's Protein Prepare module, and the protonation states of hydroxyl, Asn, Gln, and His in the structure were assigned by Schrodinger's ProtAssign module [69].
PSB36 and ZM241385 were docked, and the resulting docking poses had RMSD values less than 1.8 Å from the experimental structure (S14 Fig and S4 Table). Then, we docked C8, C9 and C10 using the same methods and chose the docking pose with the lowest binding free energy. The 6 ligand-receptor complexes C8-A 1  The MD simulations were carried out in triplicate using the PMEMD algorithm of AMBER 18 software [73]. The AMBER ff14SB force field [74] was used for the A 1 AR and A 2 AR systems. The standard protonation states of the protein residues were set at appropriate protonation states by the LEAP plugin in AMBER 18 as calculated using the H++ program [75]. The GAFF force field [76] and AMBER's Antechamber were used to generate the topology and coordinate files of the compounds C8, C9 and C10. Lipid14 force field [77] was used for the POPE bilayer. First, each system was minimized by 10,000 steps. Second, a Langevin thermostat [78] was used to heat each system from 0 K to 310 K within 500 ps. Then, a 5-ns MD simulation in the NVT ensemble was performed while retaining the heavy atoms of the protein, ligand, and lipid head groups with a constraint of 50 kcal�mol −1 �Å −2 . Third, an MD simulation in the NPT ensemble was performed for 30 ns with a constraint for the protein and ligand. The restriction was gradually reduced from 50 (20 ns) to 10 (5 ns) and 2 kcal�mol −1 � Å −2 (5 ns). Finally, each system was simulated for 100 ns without constraints under constant pressure (NPT ensemble) using a Berendsen barostat [79]. The cutoff value for the nonbonded interactions was set to 12 Å. The SHAKE algorithm [80] was used to constrain the covalent bonds involving hydrogen atoms, and the PME algorithm [81] was applied to address remote electrostatic interactions. In the process of the dynamics simulation, the time steps were set to 2 fs, and the frames were saved for analysis every 5,000 steps. The CPPTRAJ tool in AMBER 18 and VMD software [71] were used to analyze the trajectories.
The binding free energy (ΔG eff ) [82,83] between the ligand and the protein was calculated using the MM-GBSA method [84][85][86] of the Python script MMPBSA.py [87]. For the calculation, the dielectric constant of the solvent was set to 80, and the dielectric constant of the solute (E in ) was set to 1 for the lipophilic binding area. The polar part of the free energy of desolvation (ΔG GB ) was calculated using a modified GB model developed by Onufriev et al. [88]. The nonpolar part of the desolvation free energy (ΔG SA ) was calculated based on the solvent accessible surface (SASA) prediction calculated by the LCPO algorithm [89].