Agonists

a) Identification of differentially expressed signaling genes. To infer the cell state-specific communications, we first identified differentially expressed signaling genes across all cell groups within a given scRNA-seq dataset, using the Wilcoxon rank sum test with the significance level of 0.05.
b) Calculation of ensemble average expression. To account for the noise effects, we calculated the ensemble average expression of signaling genes in a given cell group using a statistically robust mean method: $\mathrm{EM}=\frac{1}{2}{Q}_2+\frac{1}{4}\left(Q_1+Q_3\right)$ where Q₁, Q₂, and Q₃ is the first, second and third quartile of the expression levels of a signaling gene in a cell group.
c) Calculation of intercellular communication probability. We modeled ligand-receptor mediated signaling interactions using the law of mass action. Since the physical process of ligand-receptor binding involves protein-protein interactions, we used a random walk based network propagation technique^{83 (link),84 (link)} to project the gene expression profiles onto a high-confidence experimentally validated protein-protein network from STRINGdb^{83 (link),85 (link)}. Based on the projected ligand and receptor profiles, the communication probability P_i,j from cell groups i to j for a particular ligand-receptor pair k was modeled by: $P_{i,j}^k=\frac{L_i{R}_j}{K_h+L_i{R}_j}\times \left(1+\frac{A{G}_i}{K_h+A{G}_i}\right)\cdot \left(1+\frac{A{G}_j}{K_h+A{G}_j}\right)\times \frac{K_h}{K_h+A{N}_i}\cdot \frac{K_h}{K_h+A{N}_j}\times \frac{n_i{n}_j}{n^2},L_i=\sqrt[m1]{L_{i,1}\cdots L_{i,m1}},R_j=\sqrt[m2]{R_{j,1}\cdots R_{j,m2}}\cdot \frac{1+R{A}_j}{1+R{I}_j}.$ Here L_i and R_j represent the expression level of ligand L and receptor R in cell group i and cell group j, respectively. The expression level of ligand L with m1 subunits (i.e., $L_{i,1},\cdots ,L_{i,m1}$ ) was approximated by their geometric mean, implying that the zero expression of any subunit leads to an inactive ligand. Similarly, we computed the expression level of receptor R with m2 subunits. In addition, co-stimulatory and co-inhibitory membrane-bound receptors are capable of modulating signaling via the control of receptor activation^{86 (link)}. For the ligand-receptor pair with multiple co-stimulatory receptors, we computed the average expression of these co-stimulatory receptors (denoted by RA) and then used a linear function to model the positive modulation of the receptor expression. For each ligand-receptor pair with multiple co-inhibitory receptors, we modeled them using the same approach. A Hill function was used to model the interactions between L and R with a parameter K_h whose default value was set to be 0.5 as the input data has a normalized range from 0 to 1. The extracellular agonists and antagonists from both sender and receiver cells are able to directly or indirectly modulate the ligand-receptor interaction^{86 (link)}. For the ligand-receptor pair with multiple soluble agonists, we computed the average expression of these agonists (denoted by AG) and then used a Hill function to model the positive modulation of the ligand-receptor interaction. For the ligand-receptor pair with multiple soluble antagonists, we modeled them using the same approach. The effect of cell proportion in each cell group was also included in the probability calculation when analyzing unsorted single-cell transcriptomes, where n_i and n_j are the numbers of cells in cell groups i and j, respectively, and n is the total number of cells in a given dataset. Together, the communication probabilities among all pairs of cell groups across all pairs of ligand-receptor were represented by a three-dimensional array P (K × K × N), where K is the number of cell groups and N is the number of ligand-receptor pairs or signaling pathways. The communication probability of a signaling pathway was computed by summarizing the probabilities of its associated ligand-receptor pairs. It should be noted that we did not perform normalization along the second dimension of P such that $\sum _j{P}_{i,j}^k=1$  because the normalized data are not suitable for comparing the communication probability between different cell groups across multiple signaling pathways. The communication probability here only represents the interaction strength and is not exactly a probability.
d) Identification of statistically significant intercellular communications. The significant interactions between two cell groups are identified using a permutation test by randomly permuting the group labels of cells, and then recalculating the communication probability P_i,j between cell group i and cell group j through a pair of ligand L and receptor R. The p-value of each P_i,j is computed by: $p-value=\frac{\left\{\#m\mathrm{\mid }{P}_{i,j}^{\left(m\right)}\le P_{i,j},m=1,2,\cdots ,M\right\}}{M}$ where the probability P_i,j^(m) is the communication probability for the m-th permutation. M is the total number of permutations (M = 100 by default). The interactions with p-value <0.05 are considered significant.

Compound collection. The current Tox21 compound collection consists of 2,870 compounds: 1,408 provided by the NTP (Xia et al. 2008 (link)) and 1,462 provided by the U.S. EPA (Huang et al. 2009 (link); Judson et al. 2009 (link)). The structures and annotations of these compounds are publicly available (Huang 2010 ; PubChem 2007 , 2009 ). The compounds were dissolved in dimethyl sulfoxide (DMSO) and plated in 1,536-well plate format at 14 or 15 concentrations ranging from 0.1 μM to 20 mM. See Supplemental Material for more details (doi:10.1289/ehp.1002952).
β-Lactamase reporter gene assay and qHTS. GeneBLAzer β-lactamase (bla) HEK 293T cell lines that constitutively co-express a fusion protein composed of the LBDs of related human NRs coupled to the DNA-binding domain of the yeast transcription factor GAL and cell culture reagents were obtained from Invitrogen (Carlsbad, CA). See Supplemental Material for more details (doi:10.1289/ehp.1002952).When activated, these fusion proteins then stimulate bla reporter gene expression. Compound formatting and qHTS were performed as described previously (Xia et al. 2009b (link)). Briefly, the bla cells with different NRs were dispensed in 1,536-well plates for screening. After cells were incubated for 5–6 hr, compounds at 14 or 15 concentrations from the NTP and U.S. EPA collections were transferred to the assay plate with the final concentrations ranging from 0.5 nM to 92 μM. The plates were incubated for 16–18 hr at 37°C before detection mix was added, and the plates were then incubated again at room temperature for 1.5–2 hr. Fluorescence intensity (405 nm excitation, 460- and 530-nm emission) was measured using an Envision plate reader (PerkinElmer, Shelton, CT). Data were expressed as the ratio of 460-nm to 530-nm emissions. The assay performance was assessed by plate statistics (signal-to-background ratio, Z´-factor, coefficient of variation) (Zhang et al. 1999 (link)) (see Supplemental Material, Table 2).
qHTS data analysis. Analysis of compound concentration–response data was performed as previously described (Inglese et al. 2006 (link)). See Supplemental Material for more details (doi:10.1289/ehp.1002952). Briefly, raw plate reads for each titration point were first normalized relative to the positive control compound and DMSO-only wells and then corrected by applying an NCGC in-house pattern correction algorithm using compound-free control plates (i.e., DMSO-only plates) at the beginning and end of the compound plate stack. Concentration–response titration points for each compound were fitted to a four-parameter Hill equation (Hill 1910 ), yielding concentrations of half-maximal activity (AC₅₀) and maximal response (efficacy) values. Compounds were designated as class 1–4 according to the type of concentration–response curve observed (Inglese et al. 2006 (link)). Curve classes are heuristic measures of data confidence, classifying concentration–responses on the basis of efficacy, the number of data points observed above background activity, and the quality of fit. To facilitate analysis, each curve class was combined with an efficacy cutoff and converted to a numeric curve rank such that more potent and efficacious compounds with higher quality curves were assigned a higher rank (see Supplemental Material, Table 4). Curve ranks should be viewed as a numeric measure of compound activity. Compounds that showed activation/inhibition in both the ratio and the 460-nm readouts were defined as activators/inhibitors. Among the activators/inhibitors, compounds with curve rank ≥ 5 or ≤ –5 in the ratio readout were further defined as active activators (agonists)/inhibitors (antagonists). Compounds with curve rank 0 (curve class 4) in both the ratio and 460-nm readouts were defined as inactive, and compounds with other phenotypes were defined as inconclusive. Curve rank, potency, and efficacy data generated on all compounds and assays can be downloaded from Huang (2010) .

The following PRR agonists were used as controls in various experiments: synthetic triacylated lipopeptide CysSerLys₄ (Pam₃CSK₄, a TLR1/TLR2 agonist, InvivoGen); Lipopolysaccharide from Escherichia coli K12 strain (LPS-EK) and Salmonella minnesota (LPS-EM) (both ultrapure TLR4 agonists, InvivoGen); Lipoteichoic acid (LTA) from Staphylococcus aureus (TLR2/TLR6 agonist, Sigma-Aldrich); CL264 (TLR7 agonist, InvivoGen); zymosan from Saccharomyces cerevisiae (zymosan crude/ZymC, TLR2 and Dectin-1 agonist, Sigma-Aldrich); zymosan depleted from S. cerevisiae (zymosan purified/ZymP, InvivoGen); and laminarin from Laminaria digitata (Dectin-1 ligand, Sigma-Aldrich). The PRR agonists were used alone or in combination with 20 ng/mL mouse recombinant IFN-γ (Peprotech).

The 3D structure of the human dopamine D4 receptor was downloaded from the Protein Data Bank database [47 (link)] (PDB ID: 5WIU). It is a protein crystal in complex with the antagonist nemonapride [37 (link)]. After removing all co-crystallized small molecules, this structure served as a model for the inactive conformation of the receptor. In this work, the active conformation of the D4 receptor was used as a homology model built on a template of the active form of the D2 receptor with PDBID 6VMS. The pre-prepared one was downloaded courtesy of the GPCRdb [48 (link),49 (link),50 (link),51 (link)] just before the 2023 update, where AlphaFold2 was used [52 (link)]. The agonists and antagonists of the human D4 receptor sets of compounds were constructed using affinity results for 5997 compounds found in the ChEMBL database [53 (link),54 (link)]. From this set, 77 compounds showing agonist or partial agonist activity and 25 antagonists were extracted manually by verifying the results of the source studies. All compounds were hierarchically clustered using Molprint2D, the Tanimoto metric, and the complete cluster linking method implemented in Canvas from the Schrödinger [55 (link),56 (link)]. The centroid was selected from each cluster containing more than two members to constitute representative structures to further study. Two sets of test compounds were created, enlisting seven agonists and ten antagonists of the human dopamine D4 receptor, respectively (Figures S1 and S2).