Six distinct NFκB signaling codons convey discrete information to distinguish stimuli and enable appropriate macrophage responses.

RESOURCE AVAILABILITY

Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Alexander Hoffmann ( ahoffmann@ucla.edu ) Materials availability Mouse lines generated in this study are available upon request.

Data and code availability

All data are available at https://data.mendeley.com/datasets/6wksmvh5p4/draft?a=832656ba-2bde-40a4-8bbc-4cecb1d9543d . Software for image processing available at https://github.com/brookstaylorjr/MACKtrack . Software for computational simulations of NFκB dynamics is available at https://github.com/Adewunmi91/nfkb_model .

EXPERIMENTAL MODEL AND SUBJECT DETAILS Mouse models

The mVenus-RelA (RelA V/V ) endogenously-tagged mouse line was generated by Ingenious Targeting Laboratory. A donor sequence encoding the monomeric variant of the Venus fluorescent protein ( Koushik et al., 2006 ) joined by a short flexible linker sequence directly upstream of the start codon of the murine Rela locus was used to generate, via homologous recombination, a tagged embryonic stem cell line, that was implanted to yield heterozygous mice. These mice were then bred with a mouse line constitutively expressing the Flp recombinase to remove the Neo resistance marker included in the homologous donor sequence. We then back-crossed the resultant mice with wild-type C57BL/6J mice to remove the Flp background and generate homozygously tagged mice (RelA V/V ). mVenus-RelA mice were crossed into a IκBα −/− TNF −/+ cRel +/− line (TNF and cRel heterozygosity are required to rescue embryonic lethality of the IκBα −/− genotype) ( Shih et al., 2009 ), as well as into an IκBβ −/− IκBε −/− line ( Hoffmann et al., 2002 ). For the Sjӧgren’s syndrome mouse model, we crossed mVenus-RelA mice into a strain that harbors mutated κB sites in the IκBα promoter ( Peng et al., 2010 ).

Macrophage cell culture

Bone marrow-derived macrophages (BMDMs) were prepared by culturing bone marrow monocytes from femurs of 8–12 week old mice in CMG 14-12-conditioned medium using standard methods ( Cheng et al., 2015 ; Takeshita et al., 2000 ). BMDMs were re-plated in experimental dishes on day 4, then were stimulated on day 7. BMDMs were stimulated with indicated concentrations of lipopolysaccharide (LPS, Sigma Aldrich), murine TNF (R&D), a TLR1/2 agonist, the synthetic triacylated lipoprotein Pam3CSK4 (PAM), a TLR3 agonist, low molecular weight polyinosine-polycytidylic acid (poly(I:C) (PIC)), a TLR9 agonist, the synthetic CpG ODN 1668 (CpG).

METHOD DETAILS Biochemical assays

For immunoblots of whole cell lysates, bone-marrow derived macrophages were replated on day 4 at 20,000/cm 2 in 6-cm dishes or 6-well plates. After stimulation on day 7, sample buffer was added directly after washing cells with PBS. Immunoblots followed standard procedure with anti-RelA (sc-372, Santa Cruz Biotechnology), anti-pIKK (CST2697), and anti-IKK2 (CST2678). Western blot band intensities were quantified using ImageJ. Nuclear extract preparation and electrophoretic mobility shift assays followed published procedures ( Caldwell et al., 2014 ).

Live-cell imaging

Bone-marrow macrophages were replated on day 4 at 20,000 or 15,000/cm 2 in an 8-well ibidi SlideTek chamber, for imaging at an appropriate density (approx. 60,000/cm 2 ) on day 6 or day 7. 2 h prior to stimulation, cells were incubated for 5 min at room temperature in a solution of 2.5 ng/mL Hoechst 33342 in PBS, then BMDM culture media was replaced. This staining condition was optimized to ensure no loss of cell viability and no aberrant morphological changes over a 24 h period of imaging in the conditions described below. Cells were imaged at 5-min intervals on a Zeiss Axio Observer platform with live-cell incubation, using epifluorescent excitation from a Sutter Lambda XL light source. Images were recorded on a Hamamatsu Orca Flash 2.0 CCD camera. After the start of imaging, additional culture media containing stimulus (TNF, LPS, poly(I:C), CpG, or Pam3CSK4) was injected into the chamber in situ . We have documented the reliability of the imaging workflow by establishing that distinct biological replicates give reproducible data ( Figure S1D ) and that distinct imaging frames of the same well provide reproducible data ( Figure S1E ). All data are available at https://data.mendeley.com/datasets/6wksmvh5p4/draft?a=832656ba-2bde-40a4-8bbc-4cecb1d9543d .

Measurement of TNF secretion and surface TNF receptor expression

To measure TNF secretion, bone-marrow macrophages were replated on day 4 at 25,000/cm 2 in a 96-well format. On day 6, media was refreshed with 80 μL media containing indicated treatment (TNF, LPS, or CpG). Supernatants were collected from wells, in triplicate, at indicated time points, using procedures from the murine TNF alpha ELISA Ready-SET-Go! kit (eBioscience #88-7324-88). To optimize assay sensitivity, measurement was performed in a half-area 96-well plate (Corning #3690), and sample incubation was performed overnight at 4°C. Fluorescence measurements were performed using a standard spectrophotometer. To measure surface receptor expression, bone-marrow-derived macrophages were replated on day 4 at 20,000/cm 2 in 6-cm dishes. On day 6, media was refreshed with 3 mL media containing indicated treatment (TNF, LPS, or CpG). At indicated time point, media was rinsed out with cold PBS. Cells were incubated with fluorophore-conjugated antibodies for TNFR1, CD11b, and F4/80 (Bio-Legend #113005, eBioscience #11-0112-82, eBioscience #12-4801-82) and analyzed, in triplicate by flow cytometry. Antibody concentration and staining conditions were performed according to manufacturer recommendations. Stained cells were measured using an Accuri C6 Flow Cytometer (BD Biosystems). Fluorescence compensation and live/dead cell filtering was performed in FlowJo v10.

Measurement of single cell RNA-seq expression

BMDMs were generated from 12-week-old WT and Sjögren’s Syndrome mice, re-plated in experimental dishes on day 5 of differentiation, and stimulated on day 7 for 8 h with 100 ng/mL lipopolysaccharide (LPS, Sigma Aldrich), 10 ng/mL murine TNF (R&D), and 50 μg/mL low molecular weight polyinosine-polycytidylic acid (poly(I:C)), or media only (Untreated control). Cells were then lifted into suspension by incubating at 37 C for 5 min using Accutase, labeled with TotalSeqB hashtag antibodies (TotalSeq-B0305 – B0308 anti-mouse Hashtag Antibody) and pooled, and captured using the 10x single cell sequencing protocol. Cell viability was ensured to be > 90% at the time of capture. Libraries were prepared with the Chromium Single Cell 3′ GEM Kit, Version 3.1 Chemistry. Hashtag libraries made using the Chromium Single Cell 3′ Feature Barcode Library Kit. Samples were sequenced paired-end 2×50 on an Illumina NovaSeq 6000 instrument.

QUANTIFICATION AND STATISTICAL ANALYSIS Image analysis

Microscopy time-lapse images were exported for single-cell tracking and measurement in MATLAB R2016a. The tracking routines followed those used in earlier work ( Selimkhanov et al., 2014 ). Briefly, cells were identified using DIC images, then segmented, guided by markers from the Hoechst image. Segmented cells were linked into trajectories across successive images, then nuclear and cytoplasmic boundaries were saved and used to define measurement regions in other fluorescent channels, including mVenus-NFκB. Nuclear NFκB levels were quantified on a per-cell basis, normalized to image background levels, then were baseline-subtracted. Mitotic cells, as well as cells that drifted out of the field of view, were excluded from analysis. The toolboxes used for this analysis are available at https://github.com/brookstaylorjr/MACKtrack . Channel capacity calculation and signaling codon identification As there are ~9.3 × 10 16 seven-dimensional combinations of 918 features ( Table S3 ) and each channel capacity calculation takes ~90 s per combination, evaluating channel capacity of all combinations of features would take ~2.3 × 10 15 h (~2.7 × 10 11 years) to compute and is therefore is computationally infeasible. To narrow the search space, we utilized a feature selection approach. Since the channel capacities of individual features combine nonlinearly, there is no guarantee a high-ranking feature in low dimensional space will also be a subset of a high-ranking feature vector in high-dimensional space. Consequently, we utilized a forward feature selection approach that balances channel capacity rankings in lower dimensional space and diversity of candidates. Channel capacity calculations are performed on single dimensional features, ranked, and a subset of features above a threshold are selected to maximize diversity. As such 1D candidates are combined to form a set of 2D feature vectors. Channel capacity calculations are calculated on the 2D feature vectors, ranked and selected as in the 1D case. This iterative ranking and selection processes are repeated until additional dimensions offer no gain in channel capacity ( Table S4 ). Algorithmic detailed: We used Shannon’s information theoretic framework to correlate the stimulus condition to dynamical features extracted from temporal trajectories of NFκB activity. noise ↓ X → communication channel → Y X = stimulus condition Y = N F κ B dynamical features C ( Q )= I ( Y ; X ) I ( Y ; X ) = H diff ( Y ) − H diff ( Y ∣ X ) H diff ( X ) = ∑ i = 1 m q i H diff ( X = x i ) = − ∑ i = 1 m q i ∑ j = 1 n i 1 n i l o g 2 ( f ( X = x i ) ) H diff ( Y ) = − ∑ i = 1 m q i n i ∑ j = 1 n i log 2 ( f ( Y = y i j ) ) f ( Y = y ) = ∑ w = 1 m q w f ( Y = y ∣ X = x w ) H diff ( A ) = − ∑ j = 1 N a δ j log 2 ( f ( a j ) ) , where δ j = probability of observing a j f ( A ) = k N a V d z ( A ) ) d k V d = π 2 σ Γ ( d 2 + 1 ) H diff ( Y ∣ X ) = conditional entropy m = number of stimulus conditions n = number of cells in a condition q i = probability of observing a stimulus x i j = a single cell ′ s response k = number of neighbors used in kNN estimate of marginal distribution of Y d = vector dimension δ j = probability of observation Controlling for different sample sizes Jackknife resampling was used to control for different sample sizes by calculating channel capacity for differently-sized subsets and extrapolating to an infinite sample size. n c =24 Setting threshold t ← ( 1 2 ) [ 1 : 6 ] − ( 1 2 ) [ 1 : 6 ] t 1 ←0.3 If d > 6 then t ←[ t ,0.1*1 d −6 ] Fori = 1 … d Compute channel capacity by optimizing over marginal distribution of X For j = 1… k c j ← I ( x j ; Y ) q j ← argmax PX I ( x j ; Y ) Select a subset of feature vectors whose channel capacity values exceeds t i X *←{ x j | cj > t i }} Q *: = argmax PX I ( X *; Y ) Select a subset of feature vectors that maximizes diversity of marginal distributions Select feature vector that yields the maximum channel capacity x ^ ← { x j * ∣ c j = max ( c ) } , equivalently x ^ ← argmax x * I ( X ; Y ) Construct a set of feature vectors containing the x ^ and feature vectors whose marginal distributions, q j , are most orthogonal to q ^ ← argmax P X I ( x ^ ; Y ) x 1 o ← x ^ , q 1 o ← q ^ Form = 2… n c Q c : = { q ∣ q ∈ Q * q ^ ∉ Q 0 } q m o ← argmin Q c ‖ Q ° − Q c ‖ 2 x m o ← { x j * ∣ q j * ≡ q m o } Machine learning classification Construction of classification models We trained an ensemble of 100 decision trees using the fitcensemble function from the Statistics and Machine Learning Toolbox from MathWorks. Decision tree models are simple, highly interpretable, and can be displayed graphically ( James et al., 2013 ). Consequently, the decision process of the classifier can be easily interrogated. However, decision tree models have two key disadvantages: (1) mediocre prediction performance ( Caruana and Niculescu-Mizil, 2006 ) and (2) high variance due to overfitting ( James et al., 2013 ). Both disadvantages can be mitigated by aggregating an ensemble of decision trees. Empirical comparison of classification models shows that ensembles of decision trees outperform other classification algorithms across a variety of problem sets ( Caruana and Niculescu-Mizil, 2006 ). We used a bootstrap aggregation (bag) method for constructing the ensembles. Each tree in the ensemble is trained on a boot-strapped replica of the data—each replica is a random selection of the data with replacement. The predictions from the ensemble model are determined by a majority vote from each individual tree prediction. We trained the ensemble to learn the stimulus labels (TNF, Pam3CSK4, CpG, LPS, and poly(I:C)) from either the entire set of predictors (all 918 metrics, Table S6A ) or a subset of predictors termed “signaling codons” ( Table S6B ).

Decision tree parameters

To construct each decision tree, the software considers all possible ways to split the data into two nodes based on the values of every predictor. Then, it chooses the best splitting decision based on constraints imposed by training parameters, such as the minimum number of observations that must be present in a child node ( MinLeafSize ) and a predictor selection criterion. The software recursively splits each child node until a stopping criterion is reached. The stopping criteria include (1) obtaining a pure node that contains only observations from a single class, (2) reaching the minimum number of observations for a parent node ( MinParentSize ), (3) reaching a split that would produce a child node with fewer observations than MinLeafSize , and (4) reaching the maximum number of splits ( MaxNumSplits ). We used default values for MinLeafSize, MinParentSize , and MaxNumSplits : 1, 10, sample size − 1, respectively ( MathWorks, 2017 ). Loadings for classification models are listed in Table S6 . Since the standard prediction selection process at each node may be biased, we used a predictor selection technique, interaction-curvature test, which minimizes predictor selection bias, enhances interpretation of the model, and facilitates inference of predictor importance. The interaction-curvature technique selects a predictor to split at each node based on the p-value s of curvature and interaction tests. Whereas the curvature test examines the null hypothesis that the predictor and response variables are unassociated, the interaction test examines the null hypothesis that a pair of predictor variables and the response variable are unassociated. A node with no tests that yield p-value s ≤ 0.05 is not split. At each node, the predictor or pair of predictors that yield the minimum significant p-value (0.05) is chosen for splitting. To split the node, the software chooses the splitting rule that maximizes the impurity gain—difference in the impurity of the node (calculated using Gini’s diversity index) and the impurity of its children nodes ( MathWorks, 2017 ).

Evaluation

We evaluated the performance of the classifiers using 5-fold cross-validation, or out-of-bag (OoB) validation, or an independent testing dataset. The OoB validation is virtually identical to K-fold cross-validation ( Hastie et al., 2001 ) and imposes minimal computation costs. K-fold cross-validation increases the computational time by K fold. OoB is defined for bagged ensembles of decision trees ( Hastie et al., 2001 ); whereas K-fold cross-validation can be used agnostic of the classification algorithm and is ubiquitous. We used OoB validation primarily to evaluate dose prediction models, which can be computationally impractical when the number features get large and K-fold cross-validation is used. We used the following performance metrics: true positive rate (recall), positive predictive value (precision), area under the Receiver Operating Characteristic (ROC) curve, F1 score, Matthews correlation coefficient, markedness, informedness and mean classification margin ( Akosa, 2017 ; Powers, 2007 ; Vihinen, 2012 ).

Dose binary classification

A series of bagged decision trees were trained to classify no treatment controls and each stimulus (each dose of each ligand). The following hyperparameters were optimized using fitcensemble function in MATLAB: NumLearningCycles , MinLeafSize , MinParentSize , and MaxNumSplits were 33, 5, 2, and 100 respectively. The models were evaluated using 5-fold cross-validation. The performance metrics for the doses of each ligand were fitted to a polynomial curve using the fit function and poly3 parameter. Feature randomization Features were selected at random to match the number of component features in codewords feature set (11) using the randsample function in MATLAB. The regenerator used was ml fg6331_64 . The features were sampled 5 times. The performance values were averaged using arithmetic mean. Feature autoencoding We used a stacked autoencoder design with two autoencoders applied sequentially using trainAutoencoder and encode functions in MATLAB. The parameters for the first autoencoder are as follows: MaxEpoch , 400; L2WeightRegularization , 0.004; SparsityRegularization , 4; SparsityProportion , 0.15; ScaleData , false. The parameters for the second autoencoder are as follows: MaxEpoch , 100; L2WeightRegularization , 0.002; SparsityRegularization , 4; SparsityProportion , 0.1; ScaleData , false.

Analysis of single cell RNA-seq data

Reads were aligned to mm10 using the 10x Cell Ranger software, version 4.0. Data was processed using Cell Ranger count to obtain a counts matrix. Data was filtered by removing cells with fewer than 1500 features. TotalSeqB hashtag labels were assigned to cells when > 75% of the cell’s hashtag reads came from one barcode. The Seurat R package ( Stuart et al., 2019 ) was used to normalize the counts. PCA was run on scaled data, and Uniform Manifold Approximation and Projection (UMAP) was run through the Seurat R package using the top 20 principal components on WT and SS cells together. To determine which genes had high stimulus-specificity, ANOVA was performed for each gene for only the three stimulus conditions in WT and SS. Estimation of maximum mutual information was performed using the R package SLEMI ( Jetka et al., 2019 ). Machine learning was performed by training a random forest classifier, as implemented in the package CARET ( Kuhn, 2008 ), on 70% of the WT data for the three stimulus conditions, using 10-fold cross-validation repeated three times, and with the mtry parameter set to sqrt(# of features). The metric used to evaluate the trained model was Accuracy, since the classes were relatively balanced. Differentially expressed genes displayed in heatmaps were found using Wilcoxon Mann Whitney U tests on each stimulus condition versus others, and the top 20 genes from each condition were merged for display. GSEA was run using the package fastGSEA ( Korotkevich et al., 2019 ) on a list of genes ranked by the WT-SS difference in ANOVA F statistic, and motif analysis on the top 1000 ranked genes was done using HOMER ( Heinz et al., 2010 ) against a whole genome background ( Heinz et al., 2010 ; Korotkevich et al., 2019 ).

Mathematical modeling

Model structure Several related models of NFκB activation in response to TNF have been established and iteratively parameterized ( Ashall et al., 2009 ; Hoffmann et al., 2002 ; Tay et al., 2010 ), and used as a basis for modeling the NFκB response to LPS and other stimuli in immortalized cell lines with exogenously introduced (and overexpressed) fluorescent RelA ( Cheng et al., 2015 ; Kellogg and Tay, 2015 ). The model presented here to account for NFκB dynamics in primary macrophages is closely based on these previous studies, inheriting identical model topologies where possible and minimizing any changes to parameter values. Key experimental data constraints As a first step toward parameterizing our model, we quantified characteristics of oscillatory endogenous BMDM signaling. We observed only slight differences in peak periodicity and amplitude between conditions (roughly a 10-min difference in median period for the lowest dose of TNF which induced robust oscillations, 0.33 ng/mL, and the highest dose tested). We did, however, observe pronounced differences in duration as the dose of TNF is increased ( Figure S1F ). Median period was determined to generally fall within 90–95 min, in the same range of oscillations measured in other cell types ( Ashall et al., 2009 ; Tay et al., 2010 ). The oscillatory frequency appeared to be remarkably stable across an extremely broad range of induction levels. Indeed, the variation observed across single cells in a particular condition (or even within the same cell) is much smaller than any differences in oscillations observed between conditions. Even when other stimuli are considered, the “signature” first harmonic of the oscillatory subpopulation remains consistent. This consistency across a wide range of input conditions agrees, notably, with predictions made using simplified discrete delay model of the NFκB network ( Longo et al., 2013 ). These delays could plausibly arise from IkB mRNA (measured to be some 10–12 min) ( Mor et al., 2010 ) and protein processing. Biochemical assays indicate that the major difference between TNF and LPS-induced IKK activation is not in the maximum amplitude, but the duration of IKK induction ( Shih et al., 2009 ; Werner et al., 2005 ). TNF strongly but transiently activates IKK. Peak IKK activity is limited in duration by rapid internalization and degradation of the ligand-bound receptor ( Mosselmans et al., 1988 ; Watanabe et al., 1988 ; Werner et al., 2008 ). LPS-bound TLR4 is also rapidly internalized, but continues to strongly activate IKK from the endosome ( Zanoni et al., 2011 ). This difference is reflected in single-cell NFκB activation: while the speed of NFκB activation (roughly proportional to the peak of IKK activity) is similar between TNF and LPS, sustained high levels of IKK activity in response to LPS leads to higher peak activity ( Figures S5B and S5C ).

Model fitting

The model was first fit for TNFR signaling and TLR4 signaling, as prior work established mathematical models that recapitulate population level data ( Cheng et al., 2015 ; Werner et al., 2008 ). For the IKK-IκB-NFκB core module, model topology and parameters were confined to be near previously established values ( Table S7 ). We performed a multidimensional sweep of transport rates and found a narrow range of parameters that could account for the observed frequency invariance, with high IKK activity diminishing oscillatory behavior ( Figure S5D ). Subsequent fitting to representative NFκB trajectories (using rmsd as distance metric) allowed us to optimize other parameters, including the induced synthesis rate constant of IκBα and the activation rate constant of IKK. For the receptor-associate modules, we required the model to recapitulate rapid IKK de- and re-activation ( Behar et al., 2013 ), which allowed IKK responses to be both adaptive (in the case of TNF) and long duration (as in TLR4 responses). We employed a screen where repeated, random initialization of parameters (within an iteratively narrower range) was followed by their optimization via gradient descent (fmin function), fitting model simulations to representative NFκB trajectories. This two-stage sweep/fitting process was repeated until parameter values converged and fits to NFκB trajectory data could no longer be improved. To parameterize the TLR1/2, TLR3, and TLR9 associated signaling modules, we used prior estimates of each receptor’s abundance in monocytes/macrophages ( O’Mahony et al., 2008 ) to estimate synthesis and degradation rates. In many cases, receptor-ligand affinities were also known ( Leonard et al., 2008 ; Nakata et al., 2006 ; Rutz et al., 2004 ) and were therefore used to estimate association and dissociation of the receptor. The kinetics of each receptor’s association with a downstream adaptor (TRIF or MyD88) were taken from estimates from our TLR4 model. NFκB responses to TLR9 were observed to be more transient than to either TLR4 or TLR1/2, in agreement with previous data ( Caldwell et al., 2014 ) and the observed self-inactivation of TLR9 ( Lee et al., 2014b ). The software to run the model is available at https://github.com/Adewunmi91/nfkb_model .

Materials availability

Mouse lines generated in this study are available upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS Mouse models

The mVenus-RelA (RelA V/V ) endogenously-tagged mouse line was generated by Ingenious Targeting Laboratory. A donor sequence encoding the monomeric variant of the Venus fluorescent protein ( Koushik et al., 2006 ) joined by a short flexible linker sequence directly upstream of the start codon of the murine Rela locus was used to generate, via homologous recombination, a tagged embryonic stem cell line, that was implanted to yield heterozygous mice. These mice were then bred with a mouse line constitutively expressing the Flp recombinase to remove the Neo resistance marker included in the homologous donor sequence. We then back-crossed the resultant mice with wild-type C57BL/6J mice to remove the Flp background and generate homozygously tagged mice (RelA V/V ). mVenus-RelA mice were crossed into a IκBα −/− TNF −/+ cRel +/− line (TNF and cRel heterozygosity are required to rescue embryonic lethality of the IκBα −/− genotype) ( Shih et al., 2009 ), as well as into an IκBβ −/− IκBε −/− line ( Hoffmann et al., 2002 ). For the Sjӧgren’s syndrome mouse model, we crossed mVenus-RelA mice into a strain that harbors mutated κB sites in the IκBα promoter ( Peng et al., 2010 ).

Macrophage cell culture

Bone marrow-derived macrophages (BMDMs) were prepared by culturing bone marrow monocytes from femurs of 8–12 week old mice in CMG 14-12-conditioned medium using standard methods ( Cheng et al., 2015 ; Takeshita et al., 2000 ). BMDMs were re-plated in experimental dishes on day 4, then were stimulated on day 7. BMDMs were stimulated with indicated concentrations of lipopolysaccharide (LPS, Sigma Aldrich), murine TNF (R&D), a TLR1/2 agonist, the synthetic triacylated lipoprotein Pam3CSK4 (PAM), a TLR3 agonist, low molecular weight polyinosine-polycytidylic acid (poly(I:C) (PIC)), a TLR9 agonist, the synthetic CpG ODN 1668 (CpG).

METHOD DETAILS Biochemical assays

For immunoblots of whole cell lysates, bone-marrow derived macrophages were replated on day 4 at 20,000/cm 2 in 6-cm dishes or 6-well plates. After stimulation on day 7, sample buffer was added directly after washing cells with PBS. Immunoblots followed standard procedure with anti-RelA (sc-372, Santa Cruz Biotechnology), anti-pIKK (CST2697), and anti-IKK2 (CST2678). Western blot band intensities were quantified using ImageJ. Nuclear extract preparation and electrophoretic mobility shift assays followed published procedures ( Caldwell et al., 2014 ).

Live-cell imaging

Bone-marrow macrophages were replated on day 4 at 20,000 or 15,000/cm 2 in an 8-well ibidi SlideTek chamber, for imaging at an appropriate density (approx. 60,000/cm 2 ) on day 6 or day 7. 2 h prior to stimulation, cells were incubated for 5 min at room temperature in a solution of 2.5 ng/mL Hoechst 33342 in PBS, then BMDM culture media was replaced. This staining condition was optimized to ensure no loss of cell viability and no aberrant morphological changes over a 24 h period of imaging in the conditions described below. Cells were imaged at 5-min intervals on a Zeiss Axio Observer platform with live-cell incubation, using epifluorescent excitation from a Sutter Lambda XL light source. Images were recorded on a Hamamatsu Orca Flash 2.0 CCD camera. After the start of imaging, additional culture media containing stimulus (TNF, LPS, poly(I:C), CpG, or Pam3CSK4) was injected into the chamber in situ . We have documented the reliability of the imaging workflow by establishing that distinct biological replicates give reproducible data ( Figure S1D ) and that distinct imaging frames of the same well provide reproducible data ( Figure S1E ). All data are available at https://data.mendeley.com/datasets/6wksmvh5p4/draft?a=832656ba-2bde-40a4-8bbc-4cecb1d9543d .

Measurement of TNF secretion and surface TNF receptor expression

To measure TNF secretion, bone-marrow macrophages were replated on day 4 at 25,000/cm 2 in a 96-well format. On day 6, media was refreshed with 80 μL media containing indicated treatment (TNF, LPS, or CpG). Supernatants were collected from wells, in triplicate, at indicated time points, using procedures from the murine TNF alpha ELISA Ready-SET-Go! kit (eBioscience #88-7324-88). To optimize assay sensitivity, measurement was performed in a half-area 96-well plate (Corning #3690), and sample incubation was performed overnight at 4°C. Fluorescence measurements were performed using a standard spectrophotometer. To measure surface receptor expression, bone-marrow-derived macrophages were replated on day 4 at 20,000/cm 2 in 6-cm dishes. On day 6, media was refreshed with 3 mL media containing indicated treatment (TNF, LPS, or CpG). At indicated time point, media was rinsed out with cold PBS. Cells were incubated with fluorophore-conjugated antibodies for TNFR1, CD11b, and F4/80 (Bio-Legend #113005, eBioscience #11-0112-82, eBioscience #12-4801-82) and analyzed, in triplicate by flow cytometry. Antibody concentration and staining conditions were performed according to manufacturer recommendations. Stained cells were measured using an Accuri C6 Flow Cytometer (BD Biosystems). Fluorescence compensation and live/dead cell filtering was performed in FlowJo v10.

Measurement of single cell RNA-seq expression

BMDMs were generated from 12-week-old WT and Sjögren’s Syndrome mice, re-plated in experimental dishes on day 5 of differentiation, and stimulated on day 7 for 8 h with 100 ng/mL lipopolysaccharide (LPS, Sigma Aldrich), 10 ng/mL murine TNF (R&D), and 50 μg/mL low molecular weight polyinosine-polycytidylic acid (poly(I:C)), or media only (Untreated control). Cells were then lifted into suspension by incubating at 37 C for 5 min using Accutase, labeled with TotalSeqB hashtag antibodies (TotalSeq-B0305 – B0308 anti-mouse Hashtag Antibody) and pooled, and captured using the 10x single cell sequencing protocol. Cell viability was ensured to be > 90% at the time of capture. Libraries were prepared with the Chromium Single Cell 3′ GEM Kit, Version 3.1 Chemistry. Hashtag libraries made using the Chromium Single Cell 3′ Feature Barcode Library Kit. Samples were sequenced paired-end 2×50 on an Illumina NovaSeq 6000 instrument.

Key experimental data constraints As a first step toward parameterizing our model, we quantified characteristics of oscillatory endogenous BMDM signaling. We observed only slight differences in peak periodicity and amplitude between conditions (roughly a 10-min difference in median period for the lowest dose of TNF which induced robust oscillations, 0.33 ng/mL, and the highest dose tested). We did, however, observe pronounced differences in duration as the dose of TNF is increased ( Figure S1F ). Median period was determined to generally fall within 90–95 min, in the same range of oscillations measured in other cell types ( Ashall et al., 2009 ; Tay et al., 2010 ). The oscillatory frequency appeared to be remarkably stable across an extremely broad range of induction levels. Indeed, the variation observed across single cells in a particular condition (or even within the same cell) is much smaller than any differences in oscillations observed between conditions. Even when other stimuli are considered, the “signature” first harmonic of the oscillatory subpopulation remains consistent. This consistency across a wide range of input conditions agrees, notably, with predictions made using simplified discrete delay model of the NFκB network ( Longo et al., 2013 ). These delays could plausibly arise from IkB mRNA (measured to be some 10–12 min) ( Mor et al., 2010 ) and protein processing. Biochemical assays indicate that the major difference between TNF and LPS-induced IKK activation is not in the maximum amplitude, but the duration of IKK induction ( Shih et al., 2009 ; Werner et al., 2005 ). TNF strongly but transiently activates IKK. Peak IKK activity is limited in duration by rapid internalization and degradation of the ligand-bound receptor ( Mosselmans et al., 1988 ; Watanabe et al., 1988 ; Werner et al., 2008 ). LPS-bound TLR4 is also rapidly internalized, but continues to strongly activate IKK from the endosome ( Zanoni et al., 2011 ). This difference is reflected in single-cell NFκB activation: while the speed of NFκB activation (roughly proportional to the peak of IKK activity) is similar between TNF and LPS, sustained high levels of IKK activity in response to LPS leads to higher peak activity ( Figures S5B and S5C ).

Supplementary Material 1 2 3 4 5 6 7 8 9 10