Single-cell stimulus-response gene expression trajectories reveal the stimulus specificities of dynamic responses by single macrophages

Abstract

Macrophages induce the expression of hundreds of genes in response to immune threats. However, current technology limits our ability to capture single-cell inducible gene expression dynamics. Here, we generated high-resolution time series single-cell RNA sequencing (scRNA-seq) data from mouse macrophages responding to six stimuli, and imputed ensembles of real-time single-cell gene expression trajectories (scGETs). We found that dynamic information contained in scGETs substantially contributes to macrophage stimulus-response specificity (SRS). Dynamic information also identified correlations between immune response genes, indicating biological coordination. Furthermore, we showed that the microenvironmental context of polarizing cytokines profoundly affects scGETs, such that measuring response dynamics offered clearer discrimination of the polarization state of individual macrophage cells than single time-point measurements. Our findings highlight the important contribution of dynamic information contained in the single-cell expression responses of immune genes in characterizing the SRS and functional states of macrophages.

Keywords: context dependence; dynamics; gene expression trajectories; inducible gene expression; information theory; machine learning; macrophage; single cell; stimulus-response specificity; trajectory ensembles.

Figure 1.. Mechanistic models simulate the full trajectories of single-cell gene expression responses

(A) Schematic of single-cell signaling activity used as signal-dependent transcription factor (SDTF) inputs into the mathematical models of gene expression. Single-cell signaling activity trajectories were derived from measured averages (Cheng et al., 2017). (B) Mathematical models are used to simulate unique genes with different gene regulatory mechanisms (GRMs) and activation/degradation characteristics. Transcriptional activation and mRNA degradation rates are determined by three parameters, the dissociation constant ${\mathsf{K}}_{\mathrm{D}}$ (affinity ${\mathsf{K}}_{\text{TF}}=1∕K_D$), the Hill coefficient n, and the mRNA degradation rate ${\mathsf{k}}_{\text{deg}}$, as described by Cheng et al, 2017. (C) Example single-cell gene expression trajectories (scGETs) simulated from the models described in Fig. 1B. Five example genes are shown that have different gene regulatory mechanisms, with 50 single cell trajectories shown per gene. Genes were defined by the following parameters [genes shown have parameters taken from GRM cluster averages determined in Cheng, et al. 2017]: AP1-regulated gene (n=1, ${\mathsf{K}}_{\text{TF},\text{AP}1}=0.48$, mRNA half-life 30 mins); NFκB-regulated gene (n=1, ${\mathsf{K}}_{\text{TF},\text{NF}\kappa \mathrm{B}}=1.3$, mRNA half-life 30 mins); NFκB∣p38-regulated gene (n=1, ${\mathsf{K}}_{\text{TF},\text{NF}\kappa \mathrm{B}}=1.3$, mRNA half-life p38 dependent [Fig. S1C]); NFκB∣IRF-regulated gene (n=1, ${\mathsf{K}}_{\text{TF},\text{NF}\kappa \mathrm{B}}=1.3$, ${\mathsf{K}}_{\text{TF},\text{IRF}}=1.25$, mRNA half-life 30 mins, governed by OR gate regulatory logic); IRF-regulated gene (n=1, ${\mathsf{K}}_{\text{TF},\text{IRF}}=1.25$, mRNA half-life 30 mins). (D) Single-cell gene expression trajectories (scGETs) can be decomposed into dynamical features (AU – arbitrary units; scGETs normalized to the maximum over all trajectories). (E) Violin plots showing the single-cell distributions of dynamical feature values for the five example genes in Fig. 1C. (F) Gene-gene correlation of the Integral values across single cells are indicative of related gene regulatory mechanisms. Top: NFκB-regulated gene vs. IRF-regulated gene (same genes as Fig. 1C). Bottom: Same NFκB-regulated gene vs. NFκB∣IRF-regulated gene (same genes as Fig. 1C).

Figure 2.. An imputation method identifies single-cell gene expression trajectories (scGETs) that match time-series data distributions

(A) Schematic of a workflow to impute the ensemble of single-cell gene expression trajectories (scGETs) and quality control its performance. Key steps of the imputation workflow are depicted in blue. This workflow begins with measured scRNAseq data at several timepoints (grey), or for quality control, simulated single-cell gene expression trajectories (green). The resulting single cell gene expression trajectories (scGETs) are then decomposed in dynamic trajectory features which may be used to address three questions (orange) with higher precision than single timepoint measurements. (B) PCA on simulated data for selected timepoints, either five front-loaded timepoints (referred to as 5.selected), or five evenly-spaced timepoints, shows that selecting front-loaded timepoints produces smoother cell state trajectories. For simulated data, the same number of cells were present at each timepoint and k-means clustering was not used. Gray lines represent ground truth connections over time. (C) Heatmaps of the transition probability matrices weighted by Euclidean distance of every pair of cells from one timepoint to the next. Shown are the heatmaps for the first two timesteps from 0hrs to 0.5hrs, and 0.5hrs to 1hr. Blue boxes indicate lower Euclidean distance and thus higher probability linkages. (D) Random walks are weighted by inverse distance, so that some cell archetypes are passed through more frequently. Exponentiating the probability matrix decreases the randomness of the random walks and identifies the high probability linkages. (In applications of the method in the paper, the probability matrix was not exponentiated.) Two archetypes in particular are included at higher frequency at 0.5hr, and 1 archetype at 1hr, corresponding to the archetypes in panel C that have low pairwise Euclidean distance across adjacent timepoints. (E) Spline interpolation over the paths determined by weighted random walks. Top: PC1 (left) or PC2 (right) for a single-cell trajectory with respect to real-time. Colored points indicate the distribution of PC scores of single-cells gene expression values collected at each timepoint. Bottom: PC1 vs. PC2 plotted for a single cell (left) or hundreds of cells (right), given the distribution of single-cell gene expression values collected at each timepoint. (F) Recovery of single-cell trajectories for individual genes (scGETs) from the interpolated PC scores matrix, using the loadings matrix from the original PCA. Shown is an NFκB-regulated gene (gene defined by parameters n=1, ${\mathsf{K}}_{\text{TF},\text{NF}\kappa \mathrm{B}}=1.3$, mRNA half-life 30 mins). (G) Violin plots of the distribution of trajectory features obtained from decomposing imputed scGETs (illustrated in Fig. 1D and defined in Methods), for the NFκB-regulated gene in Fig. 2F and three other genes that only vary in mRNA half-life. (H) Based on Euclidian distance, weighted random walks are closer to ground truth paths than completely random connections, given an equivalent number and spacing of timepoints, across a range of timepoints used. Number of principal components is held constant at 50. Five selected front-loaded timepoints (5.selected) are closer to the ground truth than 5 evenly-spaced timepoints (5). (I) Dynamical features calculated from the ground truth model-simulated scGETs show similar distributions to those imputed (Fig. 2G).

Figure 3.. Imputation of scGETs in single macrophage cells responding to immune threats

(A) Schematic of time-series scRNAseq data collected for three macrophage states responding to six different immune stimuli, totaling over 72 scRNAseq samples (3 macrophage states x 6 stimuli x at least 4 timepoints). (B) Violin plots illustrate the stimulus-response distributions of single-cell Tnf vs. Cxcl10 expression for M0 macrophages responding to each of the labeled stimuli, at each of the measured timepoints. Black points represent single cells. Green points represent the median across cells. (C) Left: Heatmap of measured time-series data for Tnf. Once the cell is sampled, it is destroyed, leading to gray boxes (no measurement) across time. Right: Imputed scGETs for Tnf expression in stimulated M0 macrophages. scGETs are displayed as a heatmap scaled to max expression across all stimuli. Row annotation colors indicate the stimulus each cell received. (D) Hierarchical clustering of single cells using the Tnf expression trajectory. Row colors indicate the stimulus each cell received. scGETs scaled to max expression across all stimuli. (E) Hierarchical clustering of the response trajectories of five genes from each cell that have different gene regulatory strategies, which further distinguishes single cell responses to each stimulus. Row annotation colors indicate the stimulus each cell received. scGETs are scaled to max expression across all stimuli. (F) Single-cell gene expression response trajectory distributions, with distribution boundaries drawn at +/− 2 standard deviations, for three genes with different gene regulatory strategies. Vertical dotted lines at 0, 0.25, 1, 3, and 8 represent timepoints at which single-cell distributions were measured via scRNAseq. Left: NFκB-only gene. Center: NFκB gene also affected by p38 activity. Right: IRF gene also affected by NFκB activity. scGETs are scaled to max expression across all stimuli; only three stimulus condition distributions are shown for clarity. (G) Distribution of scGETs for three example genes, for three different stimulus conditions, when 2,10, or 100 cells are imputed. Far right: For LPS-stimulated cells, distribution boundaries drawn at +/− 2 standard deviations stabilize between 50 to 100 imputed cells.

Figure 4.. Dynamical features of imputed scGETs contain more stimulus information than single timepoint values

(A) Distributions of the four dynamical features of scGETs: Peak Amp, Max LFC, Integral, Speed. Each point represents a single cell’s Tnf, Nfkbia, or Cxcl10 dynamical feature value in response to the indicated stimuli. (B) Stimulus-information (quantified by max MI) within a gene when considering single time-points (green) vs. linked timepoints (blue) vs. dynamical features (red), for all individual genes (points). Max MI values for each gene are connected by black lines across each feature type. Three known stimulus-response genes, Tnf, Nfkbia, and Cxcl10 are highlighted as examples. (C) Paired boxplots of the max MI for a gene at 3hrs post-stimulation vs. the corresponding information gain by considering the trajectory Integral. Genes were categorized by their reported gene regulatory mechanism, resulting in five categories as also mathematically modeled in Figure 1. The trajectory Integral contains more stimulus information than the 3hr timepoint across most genes. Genes regulated by NFκB∣IRF are less reliant on the temporal dynamics of gene expression; the estimated SRS of NFκB∣IRF genes does not increase as much when considering trajectory Integrals.

Figure 5.. Dynamical features of imputed scGETs identify gene-gene correlations more reliably than single timepoint scRNAseq data

(A) Scatterplot of Tnf vs. Nfkbia expression from scRNAseq data, at three of the timepoints measured after stimulation. Colors represent the stimulus the cell received. Only weak correlations are observed between these two NFκB-regulated genes at any of the response timepoints. (B) Scatterplot of Tnf vs. Nfkbia using the dynamical feature Integral (total activity) reveals much stronger correlations. The dynamical feature Integral reveals the known regulation of Tnf by NFκB∣p38 while Nfkbia is NFκB-only. The activation of MAPKp38 by the bacterial ligands LPS, CpG, and P3C result in a shift of Tnf expression from the identity line. (C) Scatterplot of Cxcl10 vs. Cxcl9 expression from scRNAseq data, at three response timepoints. Single-cell within-stimulus and across-stimulus correlations are timepoint-dependent and responses appear uncorrelated. (D) Scatterplot of Cxcl10 vs. Cxcl9 using the Integral (total activity) values from scGETs reveals within-stimulus correlations of gene expression across single cells, and across-stimulus response-specificity. Co-regulated genes within single cells can appear uncorrelated at timepoints due to different timing of activation, but Integral reveals their correlation. (E) Left: Heatmaps of hierarchically clustered Pearson’s correlations across all genes at the 3hr timepoint. Right: Hierarchical clustering of Pearson’s correlations across all genes, using Integral values from scGETs. (F) Pearson correlation values between the Integral values of scGETs in response to all 6 stimuli, of all gene pairs, separated by GRM. Genes of the same GRM or sharing the same transcription factors are correlated across single cells, and positive and negative correlations across GRMs are also detected.

Figure 6.. scGETs reveal how the SRS of macrophages is altered by polarization

(A) Hierarchical clustering of M1(IFNγ) cells (center) and M2(IL4) cells (right) compared to M0 (left), based on scGETs after stimulation, for a subset of cytokine and feedback regulator genes. Each row represents the imputed scGETs of a cell stimulated with the ligand indicated by the color code. (B) PC loadings obtained from performing PCA on the dynamical features Max LFC (top) or Speed (bottom) for all genes across all macrophage cells. Genes with high loading values are genes with highly variable dynamical feature values across polarization states, identifying genes with loss of specificity in either Max LFC or Speed. (C) Violin plots of Ifit1 and Mx1 Max LFC distributions for M0 vs. M1(IFNγ) macrophages illustrate that loss of SRS in the fold-change of IRF genes occurs in M1(IFNγ) macrophages. (D) Violin plots of Ifit1 and Mx1 Speed distributions for M0 vs. M2 (IL4) macrophages illustrate that loss of SRS in the activation speed of IRF genes occurs in M2(IL4) macrophages. (E) SRS is preserved in polarized macrophages. Line plot of SRS (max MI) for the most informative combination of genes, when considering scGET Integrals (solid lines) vs. scRNAseq data at the 3hr timepoint (dotted lines). (Inset displays the genes selected as most informative in combination.) When considering Integral, the best 3-gene combination can account for nearly the maximum possible SRS of macrophage, ~2.58 bits, whereas the best 3-gene combination when considering a single-timepoint post-stimulation (3hrs) does not reveal as high a SRS. (F) Scatterplots of scGET Integrals for each polarization state show that the different combinations of three complementary genes do indeed provide high SRS.

Figure 7.. The scGETs of marker genes distinguish macrophage polarization states more precisely than timepoint measurements

(A) Scatterplots of scRNAseq expression values of three canonical M1(IFNγ) or M2(IL4) marker genes at steady-state (unstimulated) vs. the LPS-response Integral of scGETs dynamics after 8 hours of LPS stimulation. The response Integral of these canonical markers better distinguishes polarization states. (B) Violin plots of time-series scRNAseq of one of the canonical M1(IFNγ) macrophage marker Nos2, across polarization states. (C) F1 score of polarization state classification accuracy when considering the three marker genes for unstimulated macrophages at steady-state (0hrs), stimulus-response timepoints (1hr, 3hr, 8hr), or scGETs stimulus-response dynamics/dynamical features (Integral, Peak Amp, Max LFC, Speed). Classification accuracy is greater for dynamical features than for either expression at single-timepoints or expression at steady-state (0hrs).