Prior studies have shown that approximately 40% of estrogen receptor positive (ER+) breast cancer (BC) patients harbor immune signaling defects in their blood at diagnosis, and the presence of these defects predicts overall survival. Therefore, it is of interest to quantitatively characterize and measure signaling errors in immune signaling systems in these patients. Here we propose a novel approach combining communication theory and signal processing concepts to model ligand discrimination in immune cells in the peripheral blood. We use the model to measure the specificity of ligand discrimination in the presence of molecular noise by estimating the probability of error, which is the probability of making a wrong ligand identification. We apply our model to the JAK/STAT signaling pathway using high dimensional spectral flow cytometry measurements of transcription factors, including phosphorylated STATs and SMADs, in immune cells stimulated with several cytokines (IFNγ, IL-2, IL-6, IL-4, and IL-10) from 19 ER+ breast cancer patients and 32 healthy controls. In addition, we apply our model to 10 healthy donor samples treated with a clinically approved JAK1/2 inhibitor. Our results show reduced ligand identification accuracy and higher levels of molecular noise in BC patients as compared to healthy controls, which may indicate altered immune signaling and the potential for immune cell dysfunction in these patients. Moreover, the inhibition of JAK1/2 produces a unique pattern of signaling dysfunction, inducing increased ligand detection error rates and reduced signal-to-noise ratios for most immune cell subtypes. These results suggest a means to improve the use of signaling kinase inhibitor therapies by identifying patients with favorable ligand discrimination specificity profiles in their immune cells.

Approximately 40% of estrogen receptor-positive breast cancer patients have problems in immune cell signaling at diagnosis, which can affect survival. This study introduces a new method using ideas from communication theory to understand how immune cells recognize signals (called ligands) in the blood. By modeling how accurately cells identify these signals despite the presence of molecular "noise," we were able to measure the likelihood of errors in signal detection in immune cells. We applied this communication model to analyze blood samples collected from breast cancer patients and healthy individuals, to reveal how cells respond to various immune-stimulating molecules. The results showed that breast cancer patients had more difficulty correctly identifying signals, suggesting their immune systems may not function properly. Additionally, when healthy cells were treated with a drug that blocks certain signaling pathways (JAK1/2), they showed similar issues in correctly identifying signals. These findings could help tailor treatments by identifying patients whose immune cells are better at signal recognition and could lead to new perspectives on the causes and consequences of immune dysfunction in breast cancer.

Communication model

A growing number of studies have used information-theoretic approaches to analyze complex single-cell data and to understand reliable communication in cell populations [9–14], and to compile a compendium of responses to cytokine stimuli [15]. Shannon's information theory [16], originally developed for digital communication systems, established fundamental limits on communication and information transmission. In Shannon's formulation, the basic problem of reliable transmission of information is stated in statistical terms, using probabilistic models for information sources and communication channels. At a high level, all communication systems have an information source (input) which sends a message, encoded or mapped into a unique signal by a transmitter. The signal is then transported from the transmitter to the receiver through a channel. Regardless of the physical medium used for transmission of the information, the main feature of the channel is that the transmitted signal is corrupted in a random manner by a variety of mechanisms, such as additive noise. The receiver decodes the corrupted signal into a decoded message, which is consumed in turn by an information sink (output). Under the principles of information theory, reliable transmission of information is possible if the information rate from the source is less than the channel capacity, which is the maximal mutual information (MI) between the input and the output of the channel.

This communication model is general enough to be applied to systems other than digital communication systems, such as a cell communication system (Fig 1). An intuitive application to cellular signaling is that the source is a cell that secretes a cell-signaling protein or cytokine (the message), recognized by a receptor (the transmitter). The receptor initiates a signaling pathway (the channel), which results in the activation of a transcription factor (the receiver). The transcription factor in turn induces gene expression (the output). When cells are stimulated with different cytokines, each cell may respond differently due to molecular noise, variations in ligand-receptor affinity, and other factors.

Here we focus on the JAK/STAT intracellular signaling axis [5], with each phosphorylation event of a STAT (pSTAT) or SMAD (pSMAD) molecule modeled as a random variable with an empirical distribution. Under a Gaussian distribution assumption, the problem of signal identification is similar to the problem of signal detection in digital communications, where the noise is assumed to be an additive, zero-mean Gaussian noise [17]. The main difference is that in digital communication systems, the set of possible transmitted signals is known, while in cellular communication systems, the signals carried by cytokines are unknown. Therefore, we must first define the cytokine-induced signals from the responses determined by the phosphorylation status of the STATs. With a zero-mean Gaussian noise assumption, the signal is given by the mean of the pSTAT responses; therefore, the Gaussian assumption is conveniently apropos. Then, the basic communication model for immune signaling in the JAK/STAT signaling pathway consists of signal generation, noisy channel transmission, and signal detection. The signal generation model maps each cytokine-induced signal to a multi-dimensional vector with elements equal to the means of the transcription factor responses. This vector defines a unique pattern of response for each cytokine treatment, which we call a codeword. The noisy channel adds a zero-mean Gaussian noise to the signaled codeword, and the detector performs optimal signal detection and identification, assuming the "correct" signal is the mean of the pSTAT responses, which is unique to each sample and stimulation. We then compare characteristics of signal detection and identification in breast cancer samples to healthy donor controls to provide insights into immune dysfunction observed in breast cancer.

Unlike previous studies that applied information-theoretic concepts such as channel capacity or mutual information to determine how much information is transmitted through signaling pathways, we adopt the perspective of detection theory to determine the specificity of signaling networks. As such, we view cytokine discrimination as the following inference problem to be solved by the cell: given the number of intracellular readout molecules, determine if an extracellular ligand is present, and decide on its type. We note that this is not a determination of a "biologically correct" response, rather, a correct determination of the noisy signal as presented to the cell, which in an in vivo context would likely depend on cell extrinsic factors.

Profiling signaling response

To characterize the specificity of immune signaling responses in estrogen receptor positive (ER+) breast cancer patients, we analyzed peripheral blood mononuclear cells (PBMCs) from 51 subjects, 32 of which were healthy donors (HD), and 19 of which were newly diagnosed with ER+ breast cancer (BC) (Tables 1, A and B in S1 Text). We systematically treated the PBMCs with 5 different cytokines alone or in combination with a JAK1/2 inhibitor ruxolitinib to study signaling responses. The cytokines included IL-2, IL-4, IL-6, IL-10, and IFNγ. After stimulation for 15 minutes, cells were fixed and stained with a panel of 27 different intra- and extracellular markers, analyzed with spectral flow cytometry, and gated into distinct cell populations. Transcription factor phosphorylation was quantified by pSTAT1, pSTAT3, pSTAT4, pSTAT5, pSTAT6, and pSMAD2/3. Cell surface markers were used to identify CD4 + T cell subsets including central memory (TCM) CD4 + CD45RA-CD27 + , effector memory (TEM) CD4 + CD45RA-CD27-, and naïve T-cells CD4 + CD45RA+CD27+ (similarly for CD8 + T cells), CD20 + B cells, naïve B cells CD19 + CD27-, classical monocytes CD14 + CD16-, and CD3-CD16 + NK cells. Only cell populations with at least 200 cells were included in the analysis. The gating strategy to identify cell populations is provided in Figs A and B in S1 Text.

Modeling ligand encoding and discrimination

Given the probabilistic nature of signal detection, we formulated the input signal (ligand) identification problem as a multi-hypothesis testing problem, where the receiver decides which of signals is the input by choosing one of the following hypotheses , : baseline (no cytokine stimulation), : stimulation with IL-4, : stimulation with IL-2, : stimulation with IL-10, : stimulation with IL-6, or : stimulation with IFNγ. The signal response space is the set of all responses (pSTAT1, pSTAT3, pSTAT4, pSTAT5, pSTAT6, pSMAD2/3) to all signals (Fig 2).

In order to compute the probability of correct ligand identification, we compute the probability of choosing the hypothesis given the true hypothesis , . Then, the overall probability of error, , is given by where is the a priori probability of the hypothesis . The optimal signal detector that minimizes the probability of error is called the maximum likelihood (ML) detector [17]. Details are provided in the Materials and Methods.

In addition to the overall probability of error, we used the signal-to-noise ratio (SNR) to characterize the specificity of the signaling response. In digital communication, SNR is a measure of fidelity of signal transmission and detection by receivers. There is a direct relationship between SNR and the probability of error in signal identification. Assuming a communication system with possible discrete signals transmitted under an additive white Gaussian noise channel, the measured response to a discrete signal is given by: where is a zero-mean Gaussian random variable with variance . The SNR is then given by where if the signals are equiprobable.

Examination of ligand identification error rates in healthy donors and ER+ breast cancer patients

We applied our communication model to calculate the probability of error and the signal-to-noise ratio, to characterize the specificity of ligand identification for each sample and cell type and compared BC to HD (Fig 3). This analysis reveals how immune cell populations differ in both the signal fidelity (SNR) and probability of signal identification error between BC and HD. Naïve CD4 + T cells and classical monocytes showed the largest error rates, with classical monocytes having lower SNR as compared to CD4 + T cells. Naïve CD8 + T cells showed lower SNR and higher error rates in BC as compared to HD. In contrast, classical monocytes show the same range of SNR for BC and HD but differ by an order of magnitude in the probability of error. Other cell types, such as naïve B cells, show overlapping SNR and error rates, with some healthy donor naïve B cells showing error rates exceeding those of breast cancer patients. NK cells showed the widest range of both SNR and error rates for both HD and BC. We observed increased error rates and reduced SNR for BC as compared to HD for nearly all immune cell subtypes (Fig 3B). SNR and graphs for CD4 + , CD8 + T, and B cell subsets are provided in Figs C-E in S1 Text.

Immune signaling error profiles

In order to compare the profile of signaling error rates across all cell types in ER+ breast cancer patients, we plotted the SNR and values together on common SNR- axes (Fig 4A). When compared directly, the cells grouped into 4 distinct regions, consisting of combinations of low and high SNR and , with classical monocytes showing low signal and high error, CD4 + T cells and central and effector memory subsets showing high signal and high error, B cells (naïve and memory) showing high signal and low error, and CD8 + T cells (including TCM and TEM) showing low signal and low error. To quantitatively compare these profiles for both BC and HD samples, we performed hierarchical clustering on the negative log transformed values (Fig 4B). This analysis reveals a clear separation between BC and HD in terms of signaling specificity (columns), with a higher overall probability of signal misidentification in BC compared to HD, as well as similarity of cell types (rows), with B, CD4 + , and CD8 + cells clustering together. This analysis reveals patterns of signal detection errors across immune cell types that are similar across HD and BC samples.

In addition to the probability of error, we examined error rates for pairs of ligands, which gives rise to a confusion matrix. For many cell types, the confusion analysis resulted in small variations for most pairs of ligands, however, the analysis revealed that IL-2 and IL-4 pair-wise error rates dominate the overall probability of error in naïve CD4 + T cells. In other words, IL-2 and IL-4 are more likely to be confused one for another compared to the other cytokine pairs. The analysis also indicates a trend in elevated pair-wise error rates between IL-2 and IL-4 in CD4 + TCM cells and NK cells (S1 Figs F-K in S1 Text). No statistically significant correlations were found between signaling error rates or SNR with patient age, tumor stage, or receptor expression.

Signaling detection alterations induced by JAK inhibition

In order to investigate the impact of signaling kinase inhibition on ligand identification error rates and signal fidelity, we treated 10 of the healthy donor samples with a clinically approved JAK1/2 inhibitor, ruxolitinib. JAK1/2 inhibition increased the identification error rate and decreased the signal fidelity (increased noise) in all immune cell subtypes, resulting in a pattern of signaling dysregulation distinct from both controls and BC samples (Fig 5A). Ruxolitinib increased the cell-type specific error rate in several cell types beyond those observed in BC samples, including classical monocytes, memory and naïve B cells, and CD8 + T cells, with the notable exception of Naïve CD4 + T cells, which showed a drastic decrease in signal fidelity but no change in error rate. JAK1/2 inhibition resulted in decreased SNR for all cell types as compared to BC and control samples (Fig 5C,D).

To compare the immune cell types to each other, we plotted the ruxolitinib samples together on common SNR- axes (Fig 5B). This analysis revealed the stark differences in effect of the ruxolitinib by cell type, with classical monocytes and CD8 + T cells (including naïve, TEM, and TCM subsets) having the largest error rates and smallest SNR. In contrast, CD4 + T cells (including naïve, TEM, and TCM subsets) showed the smallest error rate and maintained a relatively large SNR.

Ethics statement

Peripheral blood samples were obtained from breast cancer patients treated at City of Hope National Medical Center in Duarte, California, in compliance with protocols approved by the Institutional Review Board with informed consent (IRB 21368 and 19186).

Human peripheral blood samples

The study cohort consisted of individuals with newly diagnosed breast cancer, all of whom were estrogen receptor-positive (ER+), progesterone receptor-positive (PR+) and HER2/neu receptor-negative (HER2-), grade Ia-IIb, with a mean age of 57.3 years at diagnosis (range 35–76 years). All samples were obtained prior to treatment. TNM tumor staging and percent positive ki67 staining were collected for each patient. A summary of patient characteristics is provided in Table A in S1 Text. Blood was collected in EDTA-treated tubes. Peripheral blood mononuclear cells were subsequently isolated using Ficoll-Paque density gradient centrifugation (Cytiva, Marlborough, MA, USA), following the manufacturer’s protocol. The isolated PBMCs were cryopreserved in a solution containing 10% dimethyl sulfoxide (DMSO) and fetal bovine serum (FBS). Age-matched healthy control samples were acquired from the City of Hope Blood Donor Center.

Cell culture

Cryopreserved PBMCs were carefully thawed and incubated overnight (16 hours) in RPMI 1640 medium, enriched with 10% fetal bovine serum and 1% penicillin-streptomycin-glutamine (PSG), under controlled conditions (37°C, 5% CO2). Cell counts and viability assessments were conducted using a hemocytometer and trypan blue exclusion method (Sigma-Aldrich). Subsequently, the cells were cultured in a 96 deep-well plate at densities ranging from 0.5 to 1 ×ばつ 10⁶ cells/ml in fresh RPMI 1640 medium (Thermo Fisher Scientific Inc., MA, USA).

Cytokine stimulation

Following a resting period, PBMCs were cultured either untreated, stimulated with cytokine alone, or in combination with 0.1 mmol/L ruxolitinib (Cayman Chemical, Ann Arbor, MI, USA), a selective Janus kinase (JAK) 1/2 inhibitor. Cytokines included IFNγ (50 ng/ml), IL-10 (50 ng/ml), IL-2 (50 ng/ml), IL-6 (50 ng/ml), or IL-4 (50 ng/ml) (PeproTech, Rocky Hill, NJ, USA) at 37°C for 15 minutes. Following stimulation, cells were fixed with 1.5% paraformaldehyde (PFA) for 10 minutes at room temperature to preserve cellular structures and signaling intermediates. Fixed cells were then washed with phosphate-buffered saline (PBS) and permeabilized using ice-cold 100% methanol. Methanol-treated cells were stored at -80°C until further analysis. Before antibody staining, the fixed and permeabilized cells were washed three times with staining buffer (PBS supplemented with 1% fetal bovine serum).

Phospho flow cytometry

Phospho flow cytometry was performed using the following antibodies: STAT4-AF647 (clone 38/p-Stat4), CD14-APC-Cy7 (clone HCD14), CD20-AF700 (clone H1), STAT6-V450 (clone 18/pStat6), PD-L1-BV510 (clone 29E.2A3), CD3-BV570 (clone UCHT1), PD1-BV605 (clone EH12.1), CD33-BV750 (clone p67.6), CD27-BV786 (clone L128), CD45RA-BUV395 (clone HI100), CD4-BUV563 (clone SK3), CD16-BUV737 (clone 3G8), CD8-BUV805 (clone SK1), STAT3-AF488 (clone 4/p-Stat3), STAT1-Percp-Cy5.5 (clone 4a), SMAD2/3-PE (clone O72-670), Foxp3-PE-CF594 (clone 259D/C7), and STAT5-PE-Cy7 (clone 47). Antibody dilutions were prepared according to the manufacturer’s instructions and optimized through preliminary experiments to achieve optimal staining. The incubation was conducted for 45 minutes at room temperature. All antibodies were sourced from BioLegend (San Diego, CA, USA) or BD Biosciences (Franklin Lakes, NJ, USA).

Data acquisition and gating strategy

Stained cells were analyzed using a Cytek Aurora flow cytometer, equipped with lasers at 355 nm, 405 nm, 488 nm, 561 nm, and 640 nm. Compensation settings were established using single-stain controls along with a negative control. Data acquisition was conducted at a rate of 1000 events per second, with between 50,000 and 100,000 events collected per sample. Gating strategies for cell population identification are provided in Figs A and B in S1 Text

Data processing

Because we observed that pSTAT or pSMAD measurements did not always follow a Gaussian distribution, with their observed distribution tending closer to a log-normal distribution, all pSTAT/pSMAD measurements were first transformed with a log-like transformation (eq. 1). Then the transformed measurements were fit to a Gaussian mixture (GM) with two components. Each component is a multivariate normal distribution, where each variate corresponds to one of the 6 signaling molecules (5 pSTATs and pSMAD2/3). The log-transformed measurements together with knowledge of the GM parameters were fed to the optimal detector in a Monte-Carlo style simulation to obtain the probability of correct signal detection, or equivalently, the average probability of error in each analyzed signaling system. Here, one signaling system refers to the signaling pathway corresponding to one HD or BC patient and one cell type.

All pSTAT/pSMAD response data were log transformed as:

where is the measured response of the kth pSTAT/pSMAD to the ith treatment, is the number of treatments (5 cytokines and untreated case), and is a parameter found by taking the minimum over all measurements belonging to a given cell type and donor category (HD or BC). A global offset was introduced to ensure the argument of the log function is greater or equal to one, since .

Measurement data was organized in comma separated value (.csv) files, one file per subject, cytokine treatment, and cell type, provided in S1 Text. Given the large amount of data, the measurements were consolidated in structure arrays, one per HD/BC subject, cell type, and cytokine treatment. Each field in the structure array is a matrix with columns (one column for each pSTAT/pSMAD) and rows, where is the number of gated cells for the ith cytokine treatment.

Gaussian mixture modeling

GM modeling assumes that the log-transformed responses have a multi-variate Gaussian mixture distribution with two components. That is, the probability density function (PDF) of is given by:

where and are multi-variate Gaussians with means and , and covariance matrices and respectively, and is the mixing parameter. The means, covariance matrices, and mixing parameters were obtained from the log-transformed measurements using the MATLAB function "fitgmdist," which implements the Expectation-Maximization algorithm.

Probability of error computation

In order to compute the probability of error, let be the probability of choosing the hypothesis given the true hypothesis . Then, the overall probability of error, , is given by where is the a priori probability of the hypothesis . The optimal signal detector that minimizes the probability of error under the Bayes strategy is one that maximizes the probability of observed data under the hypothesis . More precisely, suppose that is the multi-variate observation and are the conditional probability density functions (PDFs) of under . Moreover, assume equiprobable hypotheses, i.e., , which is a reasonable assumption in the absence of prior information on the frequencies of . Then, the optimal detector decides if is maximized. That is, it decides if This detector is called the maximum likelihood (ML) detector [17].

To evaluate the average error probability for one sample and one cell type, all the log-transformed pSTAT/pSMAD responses for all hypotheses ( matrices of size ) were fed into a simulator that implements the ML detector and accumulates the number of detection errors over all the simulated trials (). In this simulator implementation, one trial corresponds to a single cell response for each hypothesis, i.e., a row in the matrices. The simulation loops over each hypothesis and each trial , and computes the conditional probability of the hypothesis (the measured response given the hypothesis With the GM assumption, this conditional probability is given by:

where

and are the signals corresponding to the hypothesis and and are the noise covariance matrices corresponding to the hypothesis .

Then, the ML detector selects the hypothesis such that the conditional probability is maximized:

If we declare a detection error and a counter in a error matrix is incremented at the index Otherwise, a correct decision is made and the counter at the index is incremented in the error matrix. Finally, the conditional probability , the probability of choosing the hypothesis given the true hypothesis is estimated from and the overall probability of error is given by . Note that the confusion matrix is the matrix having the element indexed by ith row and jth column equal to the conditional probabilities

Signal-to-noise ratio computation

SNR is a dimensionless ratio of signal power to noise power . We extend the SNR definition to the GM mixture model as follows. We assume a communication system with possible signals which is described by the following relationship:

That is, under the hypothesis , the transmitted signal is equal to with probability and it is equal to with probability . Similarly, under the hypothesis , the noise is with probability and it is with probability , where is a zero-mean Gaussian random variable with variance and is a zero-mean Gaussian random variable with variance .

Then, for a Gaussian mixture with two components, the SNR can be expressed as:

Finally, we generalize the above expression for the multivariate signal case, where under the hypothesis , the transmitted signal is equal to with probability and it is equal to with probability . The noise corrupting the signal is a zero-mean multivariate Gaussian noise with covariance matrix , and the noise corrupting the signal is a zero-mean multivariate Gaussian noise with covariance matrix . Under these assumptions, the SNR becomes:

All communication model simulations and analyses were performed in MATLAB. Correlation and statistical analyses were performed in Prism Graphpad 10.4.0. Note that it is customary to express SNR in decibels (dB), where the relation between linear SNR and SNR in dB is given by .

S1 Text.

Table A. Breast cancer patient characteristics. Table B. Healthy donor characteristics. Fig A. Identification of PBMCs, single cells, live cells, and T, myeloid, monocyte, and B cell populations. Fig B. Identification of T cell subsets. Fig C. Signal error and SNR for CD4 + T cell subsets. Fig D. Signal error and SNR for CD8 + T cell subsets. Fig E. Signal error and SNR for B cell subsets. Fig F. Confusion matrix for CD4 + TCM cells given by . The most likely signal confusion is between IL-4 and IL-2 cytokine stimulation, followed by INF-γ and IL-2. Fig G. Confusion matrix for CD16 + NK cells given by . The most likely signal confusion is between IL-4 and IL-2 cytokine stimulation. Fig H. Confusion matrix for CD8 + T cells given by . The most likely signal confusion is between untreated and IL-10 IL-2 cytokine stimulation, followed by IL-2 and INF-γ. Fig I. Confusion matrix for naïve B cells given by . The most likely signal confusion is between IL-6 and IL-10 cytokine stimulation, followed by IL-6 and IL-2. Fig J. Confusion matrix for classical monocytes given by . The most likely signal confusion is between IL-6 and INF-γ cytokine stimulation, followed by IL-6 and IL-10. Fig K. Confusion matrix for CD20 + B cells given by . The most likely signal confusion is between IL-10 and IL-2 cytokine stimulation. Fig L. Probability of error receptor expression correlation analysis. Receptor expression was correlated with Pe for each cell type in all HD samples (n = 32, top) and BC samples (n = 19, bottom). Statistically significant positive correlations are shown in green and negative correlations are shown in red (). Fig M. Signal-to-noise receptor expression correlation analysis. Receptor expression was correlated with SNR for each cell type in all HD samples (n = 32, top) and BC samples (n = 19, bottom). Statistically significant positive correlations are shown in green and negative correlations are shown in red (). Fig N. Correlation analysis of receptor expression rSNR on classical monocytes with overall SNR in breast cancer samples. Statistically significant correlations were found in expression levels of CD130, CD119, CD122, CD210, and CD132 on classical monocytes in breast cancer samples (n = 19). CD130, CD119, and CD122 were positively correlated, and CD210 and CD132 were negatively correlated. Because overall SNR is computed by integrating all pSTAT responses from all cytokine stimulations, these correlations do not explain differences observed between healthy donors and BC samples; moreover, classical monocytes show the lowest SNR in HD and BC samples for any immune cell subtype. Fig O. Correlation analysis of PDL and PD1 expression rSNR on NK cells with Pe and overall SNR in healthy donors and breast cancer samples. Because of the clinical relevance of PDL1 and PD1 expression on immune cells, we examined the correlation of PD(L)1 expression levels in NK cells for healthy donors (n = 32) and breast cancer (n = 19) samples. Neither PDL1 nor PD1 were consistently correlated with Pe or SNR in either healthy donors or breast cancer samples. The significance of the correlation is driven by large outliers in Pe and small outliers in SNR. Fig P. Signaling profiles of ruxolitinib as compared to BC. Statistical comparisons of A. signal-to-noise ratio and B. error rates for all immune cell types for ruxolitinib (red) and breast cancer samples (blue) shown in main text Fig 5C (unpaired t-tests adjusted for multiple comparisons, ).

https://doi.org/10.1371/journal.pcbi.1013615.s001

(DOCX)