|Home | About | Journals | Submit | Contact Us | Français|
The ability to comprehensively explore the impact of bio-active molecules on human samples at the single-cell level can provide great insight for biomedical research. Mass cytometry enables quantitative single-cell analysis with deep dimensionality, but currently lacks high-throughput capability. Here we report a method termed mass-tag cellular barcoding (MCB) that increases mass cytometry throughput by sample multiplexing. 96-well format MCB was used to characterize human peripheral blood mononuclear cell (PBMC) signaling dynamics, cell-to-cell communication, the signaling variability between 8 donors, and to define the impact of 27 inhibitors on this system. For each compound, 14 phosphorylation sites were measured in 14 PBMC types, resulting in 18,816 quantified phosphorylation levels from each multiplexed sample. This high-dimensional systems-level inquiry allowed analysis across cell-type and signaling space, reclassified inhibitors, and revealed off-target effects. MCB enables high-content, high-throughput screening, with potential applications for drug discovery, pre-clinical testing, and mechanistic investigation of human disease.
High-throughput in vitro screening has accelerated the discovery of drug candidates, but paradoxically coincides with a steep decline in the approval rate for novel molecular entities1, 2. The enormous rate of attrition as drug candidates move through clinical development can be partly attributed to the disconnect between human physiology and the in vitro screening regimen (which cannot measure efficacy in heterogeneous tissues or detect off-target toxicities2–4). If the original screening regimen more closely reflected human physiology by using human samples such as PBMCs or cancer biopsies, efficacy and toxicity could potentially be identified earlier in the development process.
High-content analysis of cellular signaling networks can provide a detailed representation of cellular state5, 6; it is often presumed that additional biologically informative reads on markers of pathways would be a desirable outcome for high-throughput screening. Compounds that target certain signaling molecules can lead to successful therapeutic outcomes7, but many compounds that target known oncogenic lesions lack clinical efficacy8. As such, the in vitro targets of a drug candidate cannot be used to accurately predict efficacy in vivo due to signaling network complexity and differences between patients or cell subpopulations of the same patient6, 7, 9–12. Therefore, high-content, single cell analysis of signaling networks in human samples during drug development could provide welcome insight to the manifold effects of drugs on cellular systems.
We propose that an ideal drug screening approach should have the following features. First, it should be based on primary human samples, with systemic behavior that resembles normal physiology and the targeted disease state. Second, subpopulation-specific system-wide signaling networks and their correlation to cell and disease phenotypes should be quantified, providing a comprehensive view of the cellular state. Third, inter-cellular communication and emergent systems properties should be evaluated and fourth and finally, measurements should be performed in high-throughput. A screening approach with these features would enable the identification of compounds with in vivo efficacy against the targeted disease and low toxicity at the earliest stage of drug discovery.
Methods have been previously developed in attempts to implement these features. Parallel enzymatic or phage display assays offer exceptional in vitro selectivity profiling13–17, but do not provide in vivo data. Cellular assays based on proliferation, apoptosis, or expression of reporter proteins approximate in vivo activity18, but drug selectivity, mechanism of action, and signaling network responses cannot be determined. Gene expression analysis19, 20 and liquid chromatography coupled to tandem mass spectrometry6, 21, 22 measure thousands of parameters, but lack throughput and single-cell resolution23, 24. High-throughput microscopy offers deep characterization of single cells23–25, but the limited number of surface and signaling molecules measured simultaneously restricts the breadth of analysis.
Fluorescence-based flow cytometry (FBFC) permits measurement of up to 12 molecules on the single cell simultaneously26-28, allowing cell subpopulations and their signaling network states to be determined simultaneously29. Drug screening applications for FBFC have been implemented by hardware30, 31 or by sample multiplexing with fluorescent cell barcoding (FCB)32, 33. With these adaptations, FBFC has become a powerful tool for drug screening and pre-clinical analysis. FBFC falls short of the ideal drug screening method described above, however, because the number of simultaneously measured parameters is limited due to spectral overlap27, hampering the comprehensive analysis of signaling network states within complex cell populations.
A recent advance in flow cytometry, mass cytometry, increases the number of parameters that can be measured, reduces overlap between measurement channels, and eliminates background autofluorescence34, 35. For mass cytometry, antibodies are labeled with isotopically pure metals36 and quantified by inductively coupled plasma mass spectrometry (ICP-MS). Current labeling techniques allow for 34 parameter measurements35. The large number of parallel measurements per cell makes mass cytometry an ideal method to assay drug candidates for cellular efficacy and selectivity.
To add high sample throughput to mass cytometry and to bring it closer to the ideal screening approach outlined above, a cell-based multiplexing technique analogous to FCB has been developed, termed MCB, where each sample is encoded by a unique combination of mass-tags before multiplexing. 96-well format MCB was used to study PBMC signaling dynamics and cell-to-cell communication, measure PBMC signaling response variability between 8 donors, and to define the impact of 27 kinase inhibitors on PBMC subpopulations. The high number of simultaneously measured parameters enabled context-specific inhibitor and cell type classification by analysis within a signaling state space, as defined by the concentration of 14 signaling proteins, rather than a single molecular readout of a signaling protein and pathway. This analysis revealed that none of the compounds tested was specific for a single cell type, that inhibitor activity and selectivity were strongly context specific, and that a topology of hierarchical relationships among PBMC cell types could be recapitulated based on signaling network responses alone.
The MCB method described here is an adaptation of FCB43, 44 that improves sample throughput, reduces antibody consumption, and ensures uniformity of the antibody stain across samples. MCB uses the molecule maleimido-mono-amide-DOTA (mDOTA) (Fig. 1a). The DOTA moiety chelates rare earth metal lanthanide(III) ions with a Kd of approximately 10−16 and the maleimide moiety rapidly reacts covalently with cellular thiol groups (Fig. 1a, Supplementary Fig. 1). Using binary combinations of seven pre-loaded mDOTA-lanthanide(III) reagents, 128 (27) combinations are possible, enabling 96-well plate multiplexing (Fig. 1b). The MCB protocol is experimentally similar to the FCB protocol (Supplementary Note 1)32, 33 and after measurement, the dataset is deconvolved with Boolean gating on the mDOTA-lanthanide(III) channels (Fig. 1c).
To test the accuracy and robustness of the MCB method, two reference cell samples differing in the abundance of phosphorylation on Tyr-696 on SH2 domain-containing leukocyte protein of 76 kDa (SLP76)(Supplementary Fig. 2) were arranged in checkerboard or striped patterns on 96-well plates for MCB analysis. After 96-well multiplexing, quantification of phosphorylation on Tyr-696 and deconvolution, the patterns were accurately recovered (Fig. 1d). The resulting Z-prime values from the checkerboard and striped 96-well plates were 0.61 and 0.58, respectively, indicating that MCB is suitably robust for high-throughput drug screening applications.
To confirm that MCB allows detection of physiological signaling events in PBMCs, and to assess subpopulation-specific signaling dynamics, we used MCB (Supplementary Fig. 3) to perform a 96-sample time-course experiment from 0 to 4 hours (Fig. 2a; Supplementary File 1, Supplementary Note 2). 14 signaling nodes and 9 cell surface markers were measured over 12 stimulation conditions (Fig. 2a–b). To identify PBMC populations (Fig. 2c) and their signaling response, Spanning-tree Progression Analysis of Density-normalized Events (SPADE)35, 37 was applied (Fig. 2b–e, Supplementary Note 2).
After interferon-alpha (IFN-α) stimulation, STAT1, STAT3, and STAT5 phosphorylation was induced in most cell types38 (Fig. 2c–d, Supplementary Fig. 4); induction peaked at 15 minutes and then declined, although elevated STAT1 phosphorylation was maintained for 4 hours in B cells and natural killer (NK) cells (Fig. 2d). Unlike STAT3 and STAT5, pre-stimulation and IFN-α-induced phosphorylation levels of STAT1 varied widely among cell types; from a 2-fold induction in CD14+ monocytes, to a 5-fold induction in other cell types (Fig. 2d). In IFN-α-stimulated T cells, STAT5 phosphorylation returned to pre-stimulation levels after initial activation, but time-dependent differences in STAT5 induction were observed in T cell subtypes (Supplementary Fig. 5).
Time course analysis by MCB also allows identification of time-dependent phenomena such as feedback regulation (Supplementary Note 3, Supplementary Fig. 6–7) and inter-cellular communication. Monocytes, which express the LPS receptor toll-like receptor 4 (TLR4)39, responded first to LPS stimulation, with phosphorylation of the canonical LPS pathway members p38, ERK, and NFκB peaking at 15–30 minutes, followed by S6 phosphorylation, which peaked after 2 hours (Fig. 2e; Supplementary Fig. 8) in agreement with previously reported results40. Cells with little or no LPS receptor expression39, including B cells, T cells, and NK cells, responded to LPS at later time points, with STAT3, STAT5, and ITK phosphorylation occurring in T cells and NK cells after 2 hours and STAT1 phosphorylation in B cells after 4 hours (Fig. 2e; Supplementary Fig. 9), which is likely due to inter-cellular communication through IL6 or other factors such as TNF-α, which are known to be released by monocytes after LPS stimulation41. These results show that the MCB method can be used to identify novel, dynamic signaling events and intercellular communication on the network scale level in complex, heterogeneous cell samples.
To assess variability in signaling between donors, eight PBMC samples were interrogated by MCB (Fig. 1, Fig 3a). As in the previous experiment, 14 signaling nodes and 9 cell surface markers were measured over 12 stimulation conditions (Fig. 3a, Supplementary Fig. 6) and analyzed by SPADE37. Samples were collected 30 minutes after stimulation, as the previous time-course experiment revealed maximal signaling response at this point for most stimulus and phosphorylation site pairs.
Cell type percentages varied between the donor PBMC samples (Supplementary Table 1). Monocytes ranged from 15 % (donor 7) to 26 % (donor 3), T cells from 29 % (donor 3) to 51 % (donor 2) and B cells ranged from 4 % (donor 2) to 11 % (donor 3). A similar range of cell percentages was also visible for the cell subtypes (Supplementary Table 1). The relative expression levels of the surface markers used for immune-phenotyping were largely similar across the donor samples (e.g. CD3, Fig. 3b), except for CD33, with donors 3 and 4 showing lower expression (Fig. 3c)42.
Despite differences in cell type abundances, cell signaling in each cell population was similar across the 8 donor samples, including S6 phosphorylation after BCR/FcR-XL induction (Fig. 3d, Supplementary File 2). Systematic evaluation in signaling response similarity between donors revealed high correlation of fold-change induction for each stimulus, phosphorylation site and cell type combination between donors (Fig. 3e), ranging from 0.67 (donor 4 vs. donor 6) to 0.93 (donor 7 vs. donor 8) (Fig. 3e). Exceptions existed: contrary to all other donors, phosphorylation on STAT5 and STAT3 was hardly induced in T cells after INF-α stimulation in donor 6 (Fig. 3f, Supplementary Fig. 10), but on STAT1 (Supplementary Fig. 11). These results show that cellular signaling is largely similar between individual donor PBMC samples, even if the cell type abundances vary.
To systematically quantify PBMC type response to kinase inhibition, 96-well MCB technique was applied to PBMCs treated with an 8-step, 4-fold dose-response titration of 27 unique small molecule inhibitors, to allow calculation of the half-maximal inhibitory concentration (IC50) and percent inhibition (Fig 4a). 12 stimulation conditions were used for 30 minutes to maximize signaling space coverage (Fig. 4b). Seven channels were used for MCB, nine for cell surface marker quantification to resolve 14 cell types (Fig. 4c) and 14 to quantify protein phosphorylation sites (Fig. 4d), covering important signaling pathways in all cell types for a network-wide signaling map (Supplementary Fig. 6). A single PBMC donor sample was used for all inhibitors to allow comparability, and one inhibitor, the JAK1/JAK2 inhibitor ruxolitinib (CAS # 941678-49-5), was tested against four donors to determine inhibitor response variability.
Per inhibitor, 18,816 phosphorylation site levels were quantified (12 stimuli × 8 doses × 14 cell types × 14 phosphorylation sites), yielding 2,352 dose-response titrations (14 cell types × 14 phosphorylation sites × 12 stimuli) for a total of 63,504 dose titrations. The extracted parameters of all dose response curves, including IC50, fold-change, percent inhibition values, the corresponding confidence intervals and Z-prime scores are given in Supplementary File 3. To visualize the data, the inhibitor IC50 values and percent inhibition were systematically organized and plotted in a cell type specific manner into a two-dimensional layout guided by canonical pathways (Fig. 5, Supplementary Files 4–6).
In order to assess inhibitor selectivity, we compared the known targets of each inhibitor (Supplementary Fig. 6) to its MCB-generated inhibition fingerprint (Fig. 5; Supplementary Files 4–6). Kinase inhibitors with a wide range of targets based on their in vitro inhibition profiles such as staurosporine(CAS # 62996-74-1)13, 14, 16, 17 (Fig. 5b–c, column 22) or the receptor tyrosine kinase (RTK) inhibitors sunitinib (CAS # 557795-19-4)14, 16, 17 (Fig. 5b–c, column 24), reduced phosphorylation levels of at least one measured signaling protein in all cell types under all conditions while other inhibitors showed more selectivity for cell type and stimulus. Comprehensive analysis of signaling network response patterns allowed putative determination of inhibitor selectivity:
JAK-STAT pathway coupled cytokine receptors activate a specific set of JAKs (of the four possible: JAKs 1–3, TYK2) upon ligand binding, which in turn phosphorylate a defined subset of STAT proteins (STAT1–6)43, 44. Seven of the stimuli used (IFN-α, IFN-γ, G-CSF, GM-CSF, IL-2, IL-3, IL-12) induce previously reported JAK-STAT pairs (Supplementary Table 1), allowing the selectivity of JAK inhibitors to be assessed. These inhibitors included ruxolitinib45, clinically approved for treatment of myelofibrosis, tofacitinib (CAS # 477600-75-2), a JAK3 inhibitor in phase 3 clinical testing against rheumatoid arthritis, lestauritinib (CAS # 111358-88-4), a JAK2 and tyrosine kinase inhibitor entering a phase 1 clinical, and several research tool compounds, including JAK2 inhibitor III (CAS # 118372-34-2), JAK3 inhibitor VI (CAS # 856436-16-3) and pan-JAK inhibitor I (CAS # 457081-03-7).
Reduced phosphorylation of various STATs after GM-CSF, IL3, IL2, G-CSF, IFN-α, and IFN-γ stimulation in diverse cell types was observed after tofacinib inhibition, indicating pan-JAK inhibition (Supplementary Fig. 12a). Higher IC50 values on JAK2 dependent phosphorylation of STAT5 after IL-3 or GM-CSF stimulation compared to JAK1 and/or JAK3 dependent phosphorylation after IL-2 or IFN-α stimulation (Supplementary Fig. 12a & 12b) suggest that both JAK1 and JAK3 and to a lesser extent JAK2 or TYK2 are inhibited by tofacinib. This is in agreement with our in vitro kinase inhibition profile (Supplementary Table 3), but differs slightly to other prior in vitro data16.
Broad inhibitory effects on STAT phosphorylation after GM-CSF, IL3, IL2, G-CSF, IFN-α, and IFN-γ stimulation were also observed for ruxolitinib (Supplementary Fig. 13), lestauritinib (Supplementary Fig. 14) and the pan-JAK inhibitor I (Supplementary Fig. 15), in agreement with our in vitro inhibition profile (Supplementary Table 3). Furthermore, lestauritinib and the pan-JAK inhibitor showed sizeable effects on signaling outside the JAK pathways (Supplementary Fig. 14 & 15), indicating that these inhibitors broadly affected many signaling network nodes.
Detailed inhibition profile analysis of JAK2 inhibitor III (Supplementary Note 4, Supplementary Fig. 16a–b) and JAK3 inhibitor VI (Supplementary Note 5, Fig. 17a–b) indicated inhibition of TYK2 activity rather than JAK2 activity and Jak1 and Jak3, respectively (Supplementary Results). Comparison of the JAK2 inhibitor III MCB results with the in vitro kinase assay results were surprising (Supplementary Table 3): it did not inhibit JAK-family kinases at concentrations up to 10µM. This discrepancy between in vitro and in vivo results could be due to an allosteric mechanism of inhibition not recapitulated in vitro, or additional off-target effects. The JAK2 Inhibitor III structure suggests that it is not an ATP-competitive inhibitor, because it is bulkier than most ATP-competitive kinase inhibitors and it lacks the critical H-bond donor and acceptor pair46.
These and analysis for inhibitors of the PI3K-AKT-mTOR-p70S6K signaling pathway (Supplementary Note 6, Fig. 5b colored boxes, Supplementary Fig. 18–20) show that the detailed analysis of inhibitor-induced signaling state by MCB provides the cellular inhibitor fingerprint and target selectivity with unprecedented resolution and throughput in complex cellular mixtures.
Simultaneous quantification of signaling responses in 14 cell types in parallel allowed analysis of cell type selectivity for each inhibitor. No inhibitor showed exclusive selectivity for a single cell type, and inhibitors with broad pathway activity such as staurosporine and sunitinib displayed little to no cell type selectivity (Fig. 5b–c, columns 22 and 24). In general, the inhibition profile of HLA-DRmid monocytes differed from those of other cell types (Supplementary File 4). Inhibitors of the Src family kinases (SFKs) and receptor tyrosine kinase (RTKs) dasatinib (CAS # 302962-49-8), LCK inhibitor (CAS # 213743-31-8), and PP2 (CAS # 172889-27-9) inhibited SFK downstream signaling components in monocytes compared to other cell types, including SYK, PLCγ2, and BLNK, often independent of stimulation (Supplementary Fig. 21–22). Many inhibitors in addition to the JAK inhibitors and sunitinib affected JAK-STAT signaling in monocytes (Supplementary Fig. 21–22), independent of stimulation conditions indicating that under the conditions of our assay, SFK and JAK-STAT signaling pathways are active in monocytes, but inactive in T cells, B cells, dendritic cells, and NK cells.
The data also enabled the comparative analysis of signaling network responses to inhibition in closely related cell types. Only few compounds differ at this level, including imatinib (Supplementary Note 7, Supplementary Fig. 23), the c-Jun N-terminal kinase (JNK) inhibitor SP600125 (CAS # 129-56-6, Supplementary Note 7, Supplementary Fig. 24) and the aminoquinone antitumor antibiotic, streptonigrin (CAS # 3930-19-6). Streptonigrin induced differential signaling responses in monocyte subtypes (Supplementary Note 7, Supplementary Fig. 25a–d) and in T cells and B cell subtypes on S6, PLCγ2, SLP76/BLNK and STAT5, often independently of the stimulation (Supplementary Fig. 25c; Supplementary Fig. 25d, columns 10–14, yellow boxes). In CD8+ T cells and IgM+ B cells streptonigrin had low IC50 values and strongly inhibited phosphorylation, but only a weak inhibition was observed in CD4+ T cells and IgM− B cells (Supplementary Fig. 25c; Supplementary Fig. 25d, columns 10–14, yellow boxes) on the same sites. An exception was when cells were stimulated with PMA/ionomycin; inhibition of most signaling molecules was detected in most cell types (Supplementary Fig. 25d, row 11, grey and red boxes). Streptonigrin interferes with cell replication, transcription, and cell respiration47, but how this might lead to the observed differences is unclear. An overview of the cell type profiles of each inhibitor tested is shown in Supplementary File 4.
Overall, these results demonstrate that MCB can be used to reveal how different cell types and their underlying signaling network states are uniquely affected by given inhibitors, underscoring the need for deep signaling profiling for supporting the development of cell type-specific compounds.
Due to the large number of signaling molecules quantified per cell type under many conditions, the inhibitor-induced cellular signaling can be used to characterize the cell types and inhibitors, ranging from the analysis of cell types and inhibitors over the complete dataset to analysis of a drug’s effect in a single cell type and condition. To analyze cell type similarity, a matrix of IC50 values representing inhibitor impact for each cell type, stimulation, and phosphorylation site was generated. Likewise, to analyze inhibitor similarity, a matrix of IC50 values representing cell type response for each inhibitor, stimulation, and phosphorylation site was created. Principal Component Analysis (PCA) was performed on each of these matrices to ask which cell types or inhibitors were similar in the reduced-dimensionality space.
In the PCA over all conditions, cell types with related immune functions such as the lymphocytes were grouped closely in PC space, forming their own cluster (Fig. 6a). A second and third cluster was formed by the monocyte lineage with closely related monocyte subtypes being also closest in PC space. The corresponding bi-clustering (Supplementary Fig. 26) of all cell types recapitulated these observations, showing that at this level of signaling network resolution, the established endpoint topology of hematopoiesis39 in PBMCs was recapitulated based on signaling network response alone.
When the input matrix was restricted to the streptonigrin-induced signaling responses (Fig. 6b), a distinct picture emerged: Closely related immune cell types were differentially affected, as CD8+ T cells, IgM+ B cells, CD14− surfaceneg. and NK cells formed their own cluster in PC space, distant from CD4+ T cells and IgM− B cells. The PCA also captured the differences in inhibitor impact among the monocyte subtypes compared to the PCA over all conditions, in concordance with the previous section.
In a second analysis, the matrix was pivoted to ask which inhibitors clustered together in the reduced-dimensionality space defined by all cell type signaling states and conditions (Fig. 6c). The general kinase inhibitors staurosporine, lestauritinib and streptonigrin and the SFK and RTK inhibitors sunitinib, PP2, dasatinib and LCK inhibitor formed their own cluster in PC space, reflecting their overall high impact on signaling networks across cell types/conditions.
Restricting the input matrix to monocytes after IFN-α stimulation across all inhibitors, the question of which inhibitors similarly impact JAK-STAT signaling can be addressed (Fig. 6d). JAK2 inhibitor III, JAK3 inhibitor VI, tofacinib, crassin (CAS # 28028-68-4), SP600125 (a JNK inhibitor), and VX680 inhibitor (CAS # 639089-54-6) formed a tight group. As described above, the JAK inhibitors inhibited JAK-STAT signaling after IFN-α stimulation in most cell types, however, also Crassin displayed a similar profile (Fig. 5c, row 4, pink box) reproducing a recent finding32. SP600125, besides its assumed selectivity, inhibited phosphorylation of STAT3, STAT5, but not STAT1 (Fig. 5c, column 21, yellow box; Supplementary Fig. 24, Supplementary Fig. 27), across many cell types. The ability of SP600125 to inhibit JAK-STAT signaling was confirmed by in vitro kinase inhibition assays (Supplementary Table 3). Here IC50 values of 974nM, 736nM, 344nM and 440nM were measured for JAK1, JAK2, JAK3 and TYK2, respectively. VX680, an inhibitor of Aurora kinases, which is also active against BCR-ABL, FLT3, and JAK248, was close to SP600125 in PC component space. Inspection of Fig. 5c (column 28, blue box) showed that phosphorylation of STAT3 and STAT5 were inhibited in the presence of VX680 after IFN-α stimulation in many cell types (Supplementary Fig. 28a; Supplementary Fig. 28b, row 5, green box). This suggests that VX680 inhibits JAK1 and, even more potently according to the in vitro data, TYK2 (Supplementary Table 3). However, for VX680 no or only weak inhibition of phosphorylation on STAT5 was observed after GM-CSF or IL3 stimulation, indicating, contrary to the in vitro data, absence of JAK2 inhibition (Supplementary Fig. 28b, red boxes).
These results show that PCA allows characterization and identification of similar cell type responses to a given inhibitor and that the inhibitor-induced signaling states were sometimes independent of cell type immune function/identity, indicating enormous cellular signaling network plasticity. In addition, PCA allowed rapid classification of inhibitors based on their profiles at a given drug exposure or in a given experimental condition, and suggests novel specificities for inhibitors SP600125 and VX680.
To establish whether the inhibition datasets generated from a single PBMC donor are generalizable, or if there is variability in inhibitor response between donors, we measured the impact of ruxolitinib on four of the eight donor samples previously described that best represent the variability between donors (Fig. 3). Supplementary Fig. 29 shows that the response to inhibition between donors was similar overall, but also showed marked differences. While ruxolitinib inhibited INF-α stimulated phosphorylation on STAT1 on IgM+ B cells, IgM− B cells and CD4+ T cells in all donors analyzed (Supplementary Fig. 29, green boxes), the same site was only inhibited in two out of four donors in CD8+ T cells (Supplementary Fig. 29, row 10, red box). Similarly, G-CSF induced phosphorylation on STAT3 in CD14− HLA-DRmid monocytes was inhibited in all donors except donor 4 (Supplementary Fig. 29, row 4, blue box). Closer inspection of these differences in inhibitor response revealed they were often due to inhibition curves that fall directly above or below the R2/fold change cutoff used as a threshold for calling a site “inhibited” (see Material and Methods), and this was often compounded by differences in the level of pathway activation observed between donors after stimulation (Fig. 3). We have observed such fluctuations in human PBMCs particularly in cases of chronic diseases involving inflammation (GPN, data not shown) - indicating in part that the differences observed might indicate differing “set points” in cell subset specific activation due to prior immune encounters.
Therefore, it can be concluded that the 27 state-based kinase inhibitor profiles previously described are a comprehensive resource describing normal healthy immune response to kinase inhibition, but also underlines the need to measure several donor samples if an inhibitor must be extensively analyzed, e.g. before a clinical trial.
In the last analysis we asked how the data generated by MCB agrees with existing in vitro kinase assay data16, 17. Again we used the matrix of IC50 values representing inhibitor impact for each cell type, stimulation, and phosphorylation and the dataset from Anastassiadis et al.16 and Davis et al.17, containing kinome-wide in vitro inhibition/IC50 values for 14/9 of the compounds analyzed in this study (Supplementary File 7). For all datasets, pairwise distances between the compounds were computed. To assess the correlation between the in vivo and in vitro datasets, the pairwise distances were plotted against each other (Fig. 6e–f).
In general, inhibitors being related or unrelated to other inhibitors in the in vitro dataset displayed a similar behavior in the in vivo datasets, resulting in correlation coefficients of 0.74 and 0.75, respectively. An exception was JAK3 inhibitor VI. Here the distance was greater to most inhibitors in the in vitro dataset compared to the in vivo dataset. Nevertheless, for most compounds analyzed throughout this manuscript the in vitro inhibition profile (e.g., Go-6983, rapamycin, GDC-0941, JAK inhibitors, etc.) is largely in agreement with the MCB data (depending upon the cell type) and supports understanding the signaling network responses observed on a gross cellular level. But caution is indicated by the data insofar as lacking generality depending upon the stimulation condition, the cell types, and even the donor tested.
These results show that both approaches correlate and that the datasets are highly complementary for the mechanistic investigation of cellular signaling pathways.
MCB makes possible high-throughput experiments that were impractical to perform using FBFC or mass cytometry alone. MCB was used here to analyze PBMC signaling dynamics, cell-to-cell communication and to comprehensively profile small molecule drug regulators based on PBMC signaling network response. In these experiments, 18,816 phosphorylation levels were quantified in 1,344 cell populations from 96 multiplexed samples for each inhibitor.
By using n metal isotopes for binary-encoded MCB, 2n samples can be multiplexed. This allows over 10,000 samples to be multiplexed in a single tube with 15 channels remaining for antibody detection. At this scale, MCB becomes an attractive technique for high-throughput drug screening and genome-wide RNAi knockdown studies.
Despite several limitations (Supplementary Note 8), the approach presented here allows analysis that span from the systems-level down to single pathways and molecules. In the experiments described, high-level compound classification suggested novel molecular targets and indicated novel mechanisms of action for widely used kinase inhibitors. The ability to identify bio-active compounds such as JAK2 Inhibitor III that presumably would not have been identified by in vitro screening highlights a key advantage to the in vivo MCB approach.
As MCB enables signaling events to be monitored over time, it provides an opportunity to study signaling pathway connectivity, inhibitor impact on feedback signaling49,, and intercellular communication50. Time-resolved single cell analysis can reveal differences between immediate and subsequent indirect, adaptive effects caused by crosstalk of interwoven signaling pathways. Here, the in vivo MCB method and in vitro kinome screening methods16, 17 are complementary, because state-based analysis of activation and signaling dynamics in combination with kinase inhibitors that have been characterized in vitro is a potentially useful paradigm for the investigation of pathway mechanism, connectivity, and dynamics.
With high-throughput systems-level coverage of the most important signaling pathways and their interconnections for each cell type, the cellular signaling states induced by inhibitors could be used as a metric for pre-clinical development. Similar MCB analyses performed on defined disease samples could be used to categorize drug effects, or drug combinations, to eventually guide therapeutic strategies based on discrete knowledge of a patient’s cellular phenotypes and genotypes. Additionally, the MCB method could be used directly as a tool for personalized medicine, with the pathway activation and drug response of a patient’s in vivo or ex vivo tissue samples used to guide therapy decisions. This initiates the creation of a “minable” systems map of signaling states that reveals the modular formats that evolution uses to build a flexible immune network (thus linking signaling responses to underlying gene networks39), enabling more precise pharmacokinetic, pharmacodynamic, and complex biomarker discovery in a variety of clinical contexts.
All inhibitors and the used concentrations are shown in Table 1 of the Supplementary Material and Methods.
Two molar equivalents of maleimido-mono-amide-DOTA (Macrocyclics) were added to each metal chloride in 20 mM ammonium acetate, pH 6.0. Solutions were then immediately lyophilized and resulting solids were dissolved in DMSO to 10 mM for long-term storage at −20°C.
Human peripheral blood, collected based on an IRB-approved protocol, was obtained from the Stanford Blood Bank. The samples obtained from healthy donors were collected in heparin sulfate anticoagulant by leukapheresis and stored at room temperature for 4–6 hours. The peripheral blood mononuclear cells (PBMCs) were isolated by Ficoll-Paque density centrifugation. The isolated PBMCs were resuspended in freezing solution (90% FBS, 10% DMSO) and stored under liquid nitrogen for future use. For each use, PBMCs were thawed and then washed twice with room temperature PBMC media (RPMI, 5% FBS, 2 mM L-glutamine with 100 U/ml penicillin and 100 µg/ml streptomycin), incubated for 1 hour at 37°C in 5% CO2, and then stimulated as shown in Supplementary table 1 by the addition of IL-2, IL-3, IL-12, G-CSF, GM-CSF, interferon-α, interferon-γ, or LPS at 30 ng/ml, sodium orthovanadate at 125 µM, phorbol 12-myristate 13-acetate (PMA) at 50 nM, Ionomycin at 1 µg/ml, or a mixture of anti-IgG, anti-IgM, anti-IgK, and anti-IgL at 10 µg/ml each (BCR/FcR-XL). The stimuli used in this study are shown in Table 2 of the Supplementary Material and Methods.
Metal-labeled antibodies were prepared as described by Bendall et al.35. Briefly, antibodies were obtained in carrier-protein-free PBS and then prepared using the MaxPAR antibody conjugation kit (DVS Sciences) according to the manufacturer’s protocol. After determining the percent yield by measurement of absorbance at 280 nm, the metal-labeled antibodies were diluted in Candor PBS Antibody Stabilization solution (Candor Bioscience GmbH) for long-term storage at 4°C. Antibodies used in this study are listed in Table 3 of the Supplementary Material and Methods.
For this analysis K562 cells, a human myelogenous leukemia cell line, either untreated or treated with orthovanadate were used. Orthovanadate is a broadly active protein tyrosine phosphatase inhibitor that increases cellular tyrosine phosphorylation levels. The induction of SH2 domain-containing leukocyte protein of 76 kDa (SLP76) phosphorylation on Tyr-696 in the orthovanadate-treated cells was observed to be highly similar in multiplexed samples (Supplementary Fig. 2) compared to non-multiplexed ones, indicating that the MCB method does not alter mass cytometry measurement or introduce artifacts.
Approximately 20 million PBMCs were aliquoted into a 96-well 2-ml block. After resting for 60 minutes at 37°C, the PBMCs were stimulated with agents listed in Supplementary Table 1 for 0 minutes, 1 minutes, 5 minutes, 15 minutes, 30 minutes, 60 minutes, 120 minutes and 240 minutes.
Approximately 20 million PBMCs were aliquoted into a 96-well 2-ml block. After resting for 45 minutes at 37°C, the PBMCs were pretreated with the indicated small molecule kinase inhibitors for 15 minutes, and then stimulated with agents listed in Supplementary Table 1 for 30 minutes in the presence of the inhibitor.
At the indicated time point after stimulation, 1.6% formaldehyde (final concentration) was added to the PBMC media and cells were incubated at room temperature for ten minutes. The formaldehyde was then diluted to 0.8% with additional PBMC media, and the fixed cells were centrifuged at 600 × g for 5 minutes at 4°C. After aspirating the supernatant, the cell pellet was resuspended in ice-cold methanol and transferred immediately to −80°C for long-term storage.
PBMC samples in methanol were brought from −80°C to 4°C on ice, washed once with Cell Staining Media (CSM, PBS with 0.5% BSA, 0.02% NaN3), and then once with PBS. The cells were then resuspended in PBS, and DMSO stocks of the barcoding reagent were added at 1:100 (to 100nM final concentration). The cells were incubated at room temperature for 30 minutes, washed three times with CSM, and then pooled into a single FACS tube for staining with metal-labeled antibodies for 1 hour at room temperature. A staining volume of 300 µl was used (~1×108 cells/ml). After antibody staining, the cells were washed twice with CSM, and then incubated for 10 minutes at room temperature (or overnight at 4°C) with an iridium-containing intercalator (DVS Sciences) in PBS with 1.6% formaldehyde. The cells were then washed three times with CSM and once with PBS, diluted with water to approximately 106 cells per ml, and filtered through a 70-µm membrane just before analysis by mass cytometry.
Cells were analyzed on a CyTOF mass cytometer (DVS Sciences) at an event rate of ~500 cells per second. The settings of the instrument and the initial post-processing parameters were described previously34,35. For each barcoded sample several data files were recorded. These were concatenated using a script developed in house. The cadmium ion signals with the mass over charges of 110, 111, 112, and 114 were summed to create a single representative channel for the CD3-QDot 605 used in the mass cytometry analysis by the Flow Core package (http://www.bioconductor.org/packages/2.3/bioc/html/flowCore.html). Before gating of the cell subpopulations and determination of the IC50 values, the data were normalized as described previously35.
All analyses were performed by Reaction Biology Corporation, Malvern, Pennsylvania, U S A against active JAK1, JAK2, JAK3 and TYK2. The compounds analyzed are shown in table 2 of the main text. All kinase reactions were performed at 10µM ATP using a 10 step, 3-fold serial dilution with 10 µM as the highest compound concentration.
The cell events measured for the PBMC time-course experiment were analyzed using the software tool SPADE as described in the main text and previously35,37. All time resolved response curves for all cell types and stimuli are shown in Supplementary File 1.
The following summarizes the SPADE algorithm within the context of this time course dataset. First density-dependent down-sampling of all measured cell events to a defined target number with equalization of the representation of rare and abundant cell types was performed. The down-sampled cell events were then clustered based on expression of 9 cell surface markers (CD33, CD20, CD3, CD4, CD7, CD123, CD14, IgM, and HLA-DR) into phenotypically similar agglomerates of cells. Those agglomerates of cells phenotypically similar in 10 dimensions were connected via edges to draw a minimum spanning tree. Next, an up-sampling step was performed to assign each cell event from the initial dataset to the most representative agglomerate. Finally, the minimum spanning tree was projected in two dimensions, and the cell clusters of the tree were colored by median intensity level of a given measured parameter allowing visualization of marker expression across the entire cellular hierarchy.
The cell events measured for each inhibitor were gated according to the scheme shown in Figure 4 (main text). Each cell type was de-barcoded individually to account for differences in the distributions of barcode metals due to differing cell sizes. The de-barcoding was semi-automated for each barcode channel by creating a boundary at the minimum between the two peaks in the density estimate and then trimming 2.5% of the cells on each side of that boundary. Subsequently, each cell was sorted into its barcode well according to the 7-digit binary number assigned. The cells determined to be in the wells stimulated with BCR-XL then had their IgM gates re-drawn because BCR-XL masks the IgM epitope and shifts the IgM distribution to lower signal levels.
The dose-response curves were then computed for every combination of phosphorylation site, modulator, and cell type. This was done by fitting the arcsinh-transformed median signal value S at each dose to the sigmoidal functional form S = Top + (Bottom − Top)/(1+10^(hill*(LogIC50−LogDose))). The fits were calculated using MATLAB’s implementation of a trust region algorithm using a robust bi-square nonlinear least squares method with each point weighted by the inverse of the standard error of the mean. To determine which curves showed significant responses, the fitting scheme was first applied to five control plates of cells that were treated with DMSO but not inhibitor. Then the false positive rate was calculated for varying levels of R^2 and fold change cutoffs (Materials and Methods Fig. 1a). An individual curve was considered a responder if it exhibited a combination of R^2 and fold change that corresponded to a <1% false positive rate in the inhibitor-free plates (Materials and Methods Fig. 1b) for a given analyzed phosphorylation site (Materials and Methods Fig. 1c). All dose-response results are shown Supplementary Files 2–4 or can be viewed at www.cytobank.org/nolanlab; each curve is shown compared to the reference level and is overlaid on individual contour plots for each sample with DNA along their hidden X axes. Empty plots signify samples where there were zero cell counts.
The percent inhibitions reported for all drugs and conditions were those observed at the highest measured inhibitor concentration, regardless of whether saturation of inhibition was observed. For each curve, this was computed by dividing the difference between the fitted curve at zero dose and at the highest inhibitor concentration by the absolute value of the difference between the fitted curve at zero dose and the reference line (see Figure 1 of of the Supplementary Material and Methods).
Principal component analysis (PCA) was used to visualize the differences between various groups in the data, including all cell types, as well as the differences between all inhibitors. Features that were used for the PCA analysis consisted of all IC50 values. In addition to the overall feature matrix, PCA was run on data stratified by various subconditions, including stimulation conditions. The control replicates were combined and averaged for this analysis. Next, the pairwise Euclidean distance between all pairs of points (cell types or inhibitors) in PCA space was calculated. This distance was calculated including those principal components that recovered 90% of the total variance. K-means clustering was performed using these distances to determine subgroups of data. A best average silhouette value over 100 replicates of K-means was used to determine the potential number of clusters for each set of conditions, The final cluster number was determined via inspection of the silhouette plots. In order to simplify visualization of the overall relationships among data points, a minimum spanning tree was created for each cluster. In order to convey more information, the values of the data points mapped to PCA space are also represented. The color of the node represents each node’s location in the first principal component, and the size of the node represents the location of a node in the second principal component. In addition, a clustergram representing a biclustering for the data in each condition was generated using the same pairwise Euclidean distance matrix described previously.
IC50 values contributing the most to principal components of interest were found by determining those values most heavily weighed in each principal component, greater than 2 standard deviations from the mean. All of these methods were run using MATLAB’s existing toolkits.
All cell density plots and heat maps were created in Cytobank (www.cytobank.org, Cytobank, Inc.).
All collected data will be made publically available by publication on www.cytobank.org/nolanlab.
We would like to thank Matt Clutter, Kenneth Gibbs and Greg Behbahani for their experimental support and discussions, and Dana Pe’er and El-ad David Amir for their feedback on data analysis. B.B. was supported by fellowships of the Swiss National Science Foundation (SNF), the European Molecular Biology Organization (EMBO), and the Marie Curie IOF. E.R.Z. is supported by a fellowship from National Institute Of General Medical Sciences (F32GM093508). T.J. C. is supported by the Department of Defense (DoD) through the National Defense Science & Engineering Graduate Fellowship (NDSEG) Program, and the Stanford Graduate Fellowship in Science and Engineering. S.C.B. is supported by the Damon Runyon Cancer Research Foundation Fellowship (DRG-2017-09). G.P.N. is supported by the Rachford and Carlota A. Harris Endowed Professorship and grants from U19 AI057229, P01 CA034233, HHSN272200700038C, 1R01CA130826, CIRM DR1-01477 and RB2-01592, NCI RFA CA 09-011, NHLBI-HV-10-05(2), European Commission HEALTH.2010.1.2-1, and the Bill and Melinda Gates Foundation (GF12141-137101).