Research ArticleImmunology

Distinct single-cell signaling characteristics are conferred by the MyD88 and TRIF pathways during TLR4 activation

See allHide authors and affiliations

Science Signaling  14 Jul 2015:
Vol. 8, Issue 385, pp. ra69
DOI: 10.1126/scisignal.aaa5208


Toll-like receptors (TLRs) recognize specific pathogen–associated molecular patterns and initiate innate immune responses through signaling pathways that depend on the adaptor proteins MyD88 (myeloid differentiation marker 88) or TRIF (TIR domain–containing adaptor protein–inducing interferon-β). TLR4, in particular, uses both adaptor proteins to activate the transcription factor nuclear factor κB (NF-κB); however, the specificity and redundancy of these two pathways remain to be elucidated. We developed a mathematical model to show how each pathway encodes distinct dynamical features of NF-κB activity and makes distinct contributions to the high variability observed in single-cell measurements. The assembly of a macromolecular signaling platform around MyD88 associated with receptors at the cell surface determined the timing of initial responses to generate a reliable, digital NF-κB signal. In contrast, ligand-induced receptor internalization into endosomes produced noisy, delayed, yet sustained NF-κB signals through TRIF. With iterative mathematical model development, we predicted the molecular mechanisms by which the MyD88- and TRIF-mediated pathways provide ligand concentration–dependent signaling dynamics that transmit information about the pathogen threat.


Toll-like receptor (TLR) signaling involves a complex network of at least 12 different TLRs that engage in physical and functional interactions with various signal transduction proteins (1). Two particular adaptor proteins interact with the TLRs through a shared Toll/interleukin-1 (IL-1) receptor (TIR) homology domain, serving to transduce TLR activation downstream. This shared sequence (and function) was identified first in the adaptor myeloid differentiation marker 88 (MyD88) (2, 3) and later in the adaptor TIR domain–containing adaptor protein–inducing interferon-β (TRIF) (4). Because all TLR family members signal through pathways mediated by these two adaptors to activate specific gene expression responses, this overall structure forms a “bow tie” motif common in cellular signaling networks (1, 5). Among the TLRs, TLR4 was the first described in mammals (6), and it is the only one that uses both adaptor pathways to stimulate inflammation and innate immune responses (7).

With the help of its co-receptor MD2 (8), TLR4 recognizes bacterial lipopolysaccharide (LPS), a component of the outer membrane of Gram-negative bacteria, and activates the pleiotropic transcription factors nuclear factor κB (NF-κB) and interferon response factor 3 (IRF3). Upon binding to LPS and undergoing dimerization, TLR4 recruits MyD88 at the plasma membrane to stimulate the initial activation of the NF-κB–controlling kinase IKK [inhibitor of κB (IκB) kinase] (7). TLR4 also undergoes dynamin-dependent endocytosis and traffics to the early endosome (9), where it interacts with TRIF and its adaptor molecule (TRAM) to initiate the TRIF-dependent pathway that leads to IRF3 activation and a second wave of IKK activity (7).

The temporal profiles of NF-κB activity are stimulus-specific (10, 11) and are thought to represent a signaling code that determines downstream cellular responses (12). For example, whereas NF-κB activity is oscillatory in response to the proinflammatory cytokine tumor necrosis factor (TNF), its activity is steady in fibroblasts responding to LPS, as a result of autocrine feedback by TNF (10, 13, 14). However, in macrophages, TNF does not contribute substantially to LPS signaling (15, 16), and studies of the population average (4, 15, 17, 18) suggest that MyD88 may be responsible for an early peak in NF-κB activity (17) and that TRIF is required for a later phase (4). However, how NF-κB dynamics at the single-cell level are encoded in macrophages that respond to TLR ligands remains an open and important question.

There remains a fundamental disconnect between the robust signaling and gene expression patterns that are observed at a population level and the variability that characterizes individual cell responses (19, 20). This variability limits the capacity for reliable biochemical information transduction (21, 22), a characteristic required for mounting appropriate physiological responses to diverse external signals. We require, then, a better understanding of the origins, control, and consequences of the noise in the molecular network that determines variable NF-κB responses. Previous modeling efforts accounted for experimentally measured distributions of binary cell fates (23) or identified potential sources of noise in the reaction network relevant to NF-κB signaling, particularly in the expression of genes subjected to NF-κB–mediated feedback (2426); however, this variability has not been contextualized with noise sources in receptor-proximal signaling modules, and signaling model simulation work has not defined or used criteria in matching the measured variability in single-cell signaling (27).

Molecular mechanisms not only determine the dynamics of signaling, but also harbor potential noise sources that determine cell-to-cell variability. Broadly, variability in cellular responses arises from the following: (i) thermodynamic stochasticity in the reactions that directly control signaling, especially if small numbers of molecules are involved, and (ii) cell-to-cell differences in the abundance of network components or enzyme-catalyzed reactions naturally arise even between genetically identical cells as a result of thermodynamic stochasticity in gene expression and protein turnover processes that are extrinsic to the reactions considered within the model (28). A study of NF-κB dynamics revealed much less variability in cells that had recently divided from the same parent cell (“sibling analysis”), which suggests that extrinsic noise is likely a primary contributor to the observed variability (29); however, the sources of extrinsic noise have not yet been defined.

In the TLR4 network, TRIF signaling occurs from the endosome, which is dependent to unknown degrees on endosomal trafficking (30), whereas MyD88 signaling is initiated from the plasma membrane; thus, these signaling branches are both temporally and spatially separated (9, 31). MyD88 additionally contains a death domain (DD), which mediates homotypic interactions to form a macromolecular complex (termed the “Myddosome”) (32). These mechanisms, although identified, have not yet been integrated into a full reaction network to understand their role in TLR4 signaling. Here, we report the iterative development of a mathematical model of the TLR4–to–NF-κB signaling network in macrophages in the context of quantitative biochemical and live-cell imaging experimental studies. We reveal how specific molecular mechanisms within the MyD88 and TRIF pathways control dynamical features of the NF-κB response, as well as associated cell-to-cell variability, which together determine the capacity to transmit information about the nature and magnitude of an invading pathogen threat.


A mathematical model of TLR4 signaling predicts specific roles for MyD88 and TRIF

Within the TLR4 signaling network, three modules may be distinguished (Fig. 1A): a TLR4 module, which transduces the presence of LPS into downstream kinase activities (IKK and TBK1) through the MyD88- and TRIF-dependent pathways; an IRF module, which transduces TBK1 activity to the production of phosphorylated nuclear IRF3; and an NF-κB module, which determines the activity of NF-κB in the nucleus as a function of input IKK and IκB activities. We have previously established mathematical models for the NF-κB module (33, 34); thus, we focused here on establishing the topology, parameters, and behavior of the TLR4-proximal signaling module in macrophages, which are the primary pathogen sensors and the effectors of the innate immune response.

Fig. 1 Computational modeling of MyD88- and TRIF-dependent activation of the IKK–NF-κB pathway in response to TLR4 stimulation.

(A) The LPS-stimulated TLR4 signaling network is depicted as three regulatory modules: the TLR4 module activates the kinases IKK and TBK1 through the MyD88- and TRIF-dependent pathways; the NF-κB module (33, 34) produces nuclear NF-κB activity as a function of control of the degradation and synthesis of IκB; and the IRF3 module produces nuclear IRF activity as a function of import and export mechanisms that are controlled by TBK1. (B) BMDMs from wild-type (wt, purple), Trif−/− (red), and Myd88−/− (blue) mice were stimulated with LPS (100 ng/ml), and a time course of activation of IKK was measured. Circles, quantified IKK kinase activity associated with coimmunoprecipitates of NEMO (fig. S1A). Lines, data from model simulations. (C) BMDMs from wt mice were treated with LPS (100 ng/ml), and the relative abundance of nuclear pIRF3 protein was measured over time. Circles, quantification of the relative amounts of nuclear pIRF3 protein as assessed by Western blotting analysis (fig. S1B). Line, data from model simulations. (D) Wt BMDMs were treated with LPS (100 ng/ml) for the indicated times. Circles, quantification of the relative amounts of nuclear NF-κB as assessed by EMSA with an NF-κB–specific double-stranded oligonucleotide probe (fig. S1C). Line, data from model simulations. (E) Sample time lapse images, taken at the indicated times, of live RAW264.7 macrophages expressing EYFP-p65 treated with LPS (500 ng/ml). (F) Measured nuclear NF-κB responses to LPS of five single cells that were randomly selected from the 376 total cells tracked in the experiment shown in (E). (G) Simulation of NF-κB dynamics in wt, Trif−/− (MyD88-only), and Myd88−/− (TRIF-only) cells in response to LPS (0.1 ng/ml to 1 μg/ml, dark to bright line colors). (H) Model-predicted times of the first maximal (peak) accumulation of nuclear NF-κB with respect to LPS concentration in wt, MyD88-only, or TRIF-only cells. Dashed line indicates weak peaks (<10% of the maximum peak observed). (I) Model-predicted durations of the total activity of nuclear NF-κB ([NF-κBn] > 0.03 μM) with respect to LPS concentration in wt, MyD88-only, or TRIF-only cells. All experimental data shown in this figure are representative of three independent experiments.

With populations of bone marrow–derived macrophages (BMDMs) derived from wild-type, MyD88-deficient (Myd88−/−), and TRIF-deficient (Trif−/−) mice, we first measured IKK activation dynamics in response to different concentrations of LPS (1 and 100 ng/ml) through an established in vitro kinase assay with immunoprecipitated IKK. In Trif−/− BMDMs, IKK activity was induced early and transiently, reaching maximal activation at least 15 min earlier than that in Myd88−/− cells (Fig. 1B and fig. S1A). In contrast, the more slowly activated Myd88−/− BMDMs showed more persistent signaling, well past the hour-long window of activation that was observed in Trif−/− cells. The dynamics of IKK activity in wild-type cells showed both early activation and long duration, which are characteristics of the sum of the activities of MyD88 and TRIF. We similarly quantified the dynamics of the downstream transcription factors NF-κB, as measured by electrophoretic mobility shift assays (EMSA), and IRF3, as determined by measurement of IRF3 phosphorylation by Western blotting analysis (Fig. 1, C and D, and fig. S1, B and C).

By parameterizing ordinary differential equation (ODE) models of the TLR4 and IRF3 modules with these data (fig. S2, table S1, and Supplementary text), we were unable to obtain a good fit with the EMSA-measured NF-κB activity profiles, although the dynamics of IKK and IRF3 activation were recapitulated for all genotypes and LPS concentrations (Fig. 1, B to D, and fig. S1, D and E). Note that this model did not include autocrine cytokine mechanisms, because these play little role in macrophages responding to LPS (fig. S3) (15, 16). The key discrepancy lies in the propensity of the model to produce oscillations, which are supported by the delayed negative feedback loop mediated by IκBα (Fig. 1D and fig. S1E, bottom). That the measurements at the level of the cell population did not reveal oscillations suggested the possibility that single-cell dynamics were obscured by a high degree of cell-to-cell variability. To resolve this discrepancy, we performed live-cell imaging (Fig. 1E) of RAW264.7 cells (a mouse macrophage cell line) stably transduced with a lentivirus expressing the NF-κB subunit p65 fused to enhanced yellow fluorescent protein (EYFP-p65), which was under the control of the Rela promoter. After the cells were stimulated with a range of concentrations of LPS (500 pg/ml to 5 μg/ml), we tracked the nuclear translocation of NF-κB at the single-cell level. Responses were first evident in response to LPS (1 ng/ml) (fig. S4), and increasing concentrations of LPS resulted in oscillatory translocation patterns (Fig. 1F), in agreement with model simulations. Oscillatory dynamics in NF-κB signaling were previously reported in response to TNF (35, 36), whose signal is mediated by a different receptor module, but by the same IκB signaling mechanism that transmits signals by LPS. Thus, our results with LPS-TLR4 suggest that NF-κB oscillations are an intrinsic feature in the NF-κB–IκB signaling module.

With this model, we explored the pathway-specific roles in encoding NF-κB dynamics in response to a range of LPS concentrations (Fig. 1G, dark to bright lines). As expected, NF-κB activated solely by MyD88 (that is, in Trif−/− cells) showed no second-phase activity. Whereas the late-phase dynamics remained intact in cells expressing only TRIF, they showed a slowed and reduced first phase. The time of initial maximal activity in wild-type or Trif−/− cells was ~20 min earlier than that in Myd88−/− cells (Fig. 1H). In contrast, control of the duration of NF-κB activity was entirely TRIF-dependent, with the MyD88 pathway producing only transient NF-κB responses (Fig. 1I). Together, our results illustrated that the two TLR4-responsive pathways encoded specific aspects of NF-κB dynamics: early, transient activation of IKK by MyD88 determined the initial timing of the response, whereas slower, persistent activation of IKK by TRIF encoded a longer duration of NF-κB activity. This analysis provides a framework for studies of how the underlying molecular mechanisms in the activation of MyD88 and TRIF determine NF-κB dose responsiveness, dynamics, and cell-to-cell variability.

Signalosome formation determines MyD88-dependent signaling

MyD88 signaling is thought to involve the formation of a signaling complex, the Myddosome, which is composed of six molecules of MyD88, four molecules of IL-1 receptor–associated kinase 4 (IRAK4), and four molecules of IRAK1 (32, 37), which prompted us to examine its implicit signaling characteristics. We modeled Myddosome formation in accordance with a sequential binding proposal (32) in which the LPS-TLR4 complex (C) functions as a seed that may attract two MyD88 monomers (M), forming CM2, which in turn functions as a building block to form C2M4 and C3M6. The C3M6 complex constitutes a molecular platform that forms the macromolecular Myddosome by incorporating IRAK4 and IRAK1 (Fig. 2A and Supplementary text). As previously hypothesized (37), we found that the process was inherently cooperative, leading us to approximate the reaction scheme by Hill kinetics (fig. S5). Assuming MyD88 subunit interaction ratios (kf/kb from Fig. 2A) in a range of 0.1 to 10 resulted in a range of fitted Hill coefficients between 1.8 and 3.1 (Fig. 2B).

Fig. 2 Signalosome formation shapes the dose response of IKK and NF-κB activation.

(A) The seven-step Myddosome assembly model (see Supplementary text for details). Green numbers indicate reaction steps. kf and kb represent the forward and reverse reaction rates, respectively. (B) Fitted Hill coefficients for the different ratios of kf/kb from the Myddosome model. The red line indicates the fitted Hill coefficient value used in the full TLR4 model. (C and D) Dose responsiveness of the maximal (peak) activities of (C) IKK and (D) nuclear NF-κB as predicted by the full TLR4 model for the indicated Hill coefficients (from the Myddosome assembly model, and corresponding to different kf/kb ratios) in Trif−/− cells (gray gradient lines) and Myd88−/− cells (purple line). The control (ctrl) indicates the case in which the Hill coefficient = 1. (E) Violin plots of the amplitudes of the first peak NF-κB activity in single cells with respect to the concentration of LPS. The shapes of the plots show the distributions of total responses (normalized to the same area), whereas the points show the median responses. a.u., arbitrary unit. (F) Fraction of cells showing NF-κB activation above a given peak amplitude threshold (>0.55 a.u.) with respect to LPS concentration. (G) Simulated total IKK responses to different concentrations of LPS with Hill coefficients (n) of 1 (no cooperativity, top two panels) or 3 (bottom two panels). Left: Time courses of responses in the Trif−/− model. Right: Integrated IKK responses with respect to LPS concentration in the MyD88−/− and Trif−/− models. (H) Integrated IKK responses (area under the curve) in BMDMs from MyD88−/− and Trif−/− mice that were treated with the indicated concentrations of LPS and then were subjected to quantified kinase assays (fig. S6). The responses in (G) and (H) were normalized to total observed range. Single-cell data in (E) and (F) include all tracked cells (between 376 and 501 cells) in each condition, and each condition is representative of at least two independent experiments. Biochemical data in (H) are representative of three independent experiments.

We then used our TLR4 model to study the ramifications of the range of predicted Hill coefficients (Fig. 2, C and D). Whereas the Myd88−/− model predicted a slowly saturating peak response to increasing concentrations of LPS, maximal NF-κB activation in the Trif−/− model showed a faster switch from off to on, depending on the strength of MyD88 oligomerization. We measured the amplitude of the first peak of NF-κB activation in the single-cell responses of RAW264.7 cells (Fig. 2, E and F), and we observed an increase from minimal to almost saturated peak NF-κB activation over a range of concentrations of LPS from 1 to 10 ng/ml. These values were used to tune the effective half-maximal concentration (EC50) and Hill coefficient of the model (see Supplementary text) to provide a best-fit Hill coefficient of 3, which corresponds to a kf/kb ratio of 6.1 (Fig. 2, B to D, red line).

A second prediction of cooperative signalosome formation related to the duration of NF-κB signaling. Given a set of dose-dependent transient inputs (that is, upstream TLR activation at the plasma membrane, which was quickly turned off by receptor endocytosis), higher levels of cooperativity caused robust termination behavior that saturated not only the peak amplitude but also the total amount of active IKK over time (Fig. 2G). In a highly cooperative system (Hill coefficient of 3 in the Trif−/− model; Fig. 2G, bottom left), outputs above the activation threshold are nearly indistinguishable from one another, which results in a saturated dose-response curve for integrated IKK activity (Fig. 2G, bottom right). Integrating measured IKK activity in Myd88−/− and Trif−/− BMDMs (Fig. 2H and fig. S6) showed that the total MyD88-stimulated IKK activity saturated more quickly than did TRIF-stimulated IKK activity, which supports this prediction.

In summary, MyD88 provided for robust on-or-off control of NF-κB activity. This behavior was predicted by the mathematical model, which revealed inherent cooperativity in the pathway (Fig. 2B), and was measured by peak NF-κB translocation (Fig. 2E). Once the on state is reached, greater amounts of LPS do not further increase MyD88-mediated signaling (Fig. 2G). This dynamic behavior is partially explained by cooperative Myddosome formation, but additionally requires a robust post-induction deactivation mechanism, which may be mediated by the rapid internalization of TLR4 into endosomes, a process that limits the duration of MyD88 signaling. Thus, we turned to studying further the roles of endosome trafficking in determining the balance between the MyD88- and TRIF-dependent pathways.

Endosomal translocation and maturation shape the characteristics of both MyD88 and TRIF signaling

Ligand-bound TLR4 first induces MyD88 signaling from the plasma membrane, and then it is thought to be endocytosed to induce TRIF signaling from the early endosome (9); however, unbound TLR4 also traffics together with recycled plasma membrane to endosomes as part of normal constitutive processes (30), where it could encounter free LPS or endocytosed CD14-LPS complexes in the endosome. We sought to delineate the relative contributions of constitutive and ligand-induced endocytosis of TLR4 to NF-κB signaling (Fig. 3A).

Fig. 3 Ligand-induced endosomal transport of TLR4 mediates the transition from MyD88-dependent to TRIF-dependent signaling.

(A) The two proposed modes of generating activated TLR4 in endosomes: constitutive shuttling of receptors, which is followed by ligand binding in the endosome; and ligand-induced transport, in which endocytosis occurs after ligand binding. (B) Predicted mechanism-specific fluxes of the endosomal ligand-TLR4 complex in response to LPS (10 ng/ml and 1 μg/ml). “Constitutive” flux is the term given to the sum of LPS-TLR4 binding and unbinding in the early endosome, whereas “ligand-induced” flux is the sum of the endocytosis and exocytosis of ligand-TLR4 complexes. The net flux is the sum of both constitutive and ligand-induced flux, as well as the degradation of the ligand-bound TLR4 complex in the early endosome. (C) Model-predicted integrated flux over the short term (0 to 3 min) and the long term (0 to 4 hours) in response to the indicated concentrations of LPS. (D) Simulated nuclear NF-κB responses to the indicated concentrations of LPS (1 ng/ml to 1 μg/ml; dark to bright lines) in in silico knockout conditions in which the simulated cells can internalize TLR4 through constitutive means only (left) or through ligand-induced shuttling only (right).

We decomposed endosomal signaling by computational flux analysis in our mathematical model (see Supplementary text for details), comparing the transfer over time of the receptor-ligand complex from the plasma membrane to the endosome (activation through ligand-induced transport), to the amount of LPS binding to endosomal TLR4 (activation through constitutive shuttling). Overall, the ligand-induced transport of TLR4 proved to be dominant at all concentrations of LPS tested (Fig. 3, B and C). Indeed, apart from very high concentrations and early time points (Fig. 3C), constitutive TLR4 shuttling generally played a negative role in NF-κB activation. These model-derived conclusions were robust to a high degree of parameter perturbation (fig. S7).

We then examined—in computational simulations—whether only ligand-induced TLR4 transport, but not constitutive TLR4 shuttling, could support NF-κB activation. Unexpectedly, when either trafficking mechanism was removed from the model, the total TLR-mediated activation of NF-κB was enhanced (Fig. 3D compared to Fig. 1F). When TLR4 transport (by either mechanism) away from the plasma membrane was decreased, TRIF signaling was indeed reduced (and was almost eliminated in the constitutive shuttling-only model); however, MyD88 signaling was greatly enhanced in both cases because more TLR4 persisted at the plasma membrane at an amount greater than that required for the threshold of activation of the Myddosome (Fig. 3D). Thus, whereas ligand-induced receptor transport was essential for TRIF activation, both trafficking modalities were critical for enforcing the transience of MyD88 activation.

Whereas TLR4 transport to endosomes both limited MyD88 activity and enhanced TRIF activity, the maturation of early endosomes into late endosomes is associated with the termination of the TRIF-mediated response through degradation of TLR4 (38, 39). We measured single-phagosome maturation in RAW264.7 cells using Escherichia coli conjugated to the pH-responsive dye pHrodo Red (Fig. 4A). The fluorescence intensity of endosome-localized pHrodo increased over the pH ranges associated with the progress of endosomes through early to late stages before eventually plateauing after about 10 to 12 hours (fig. S8). We associated segmented pHrodo-containing endosome spots with single cells, and measured the increase in fluorescence in these single-cell spots over time, normalized to their final state, which was measured at 17.5 hours (Fig. 4B). We then computed the time taken for each cell to cross given specific maturation thresholds; depending on the threshold applied, a range of maturation times was generated, which resembled a normal distribution with a mean ranging from 4 to 10 hours and an SD of 1 to 2 hours (Fig. 4C). These times suggested a distribution of delays, after which endosomal ligand-TLR4 complexes underwent rapid degradation.

Fig. 4 Variable endosomal maturation shapes the termination of TRIF signaling.

(A) Example image (at t = 16 hours) of RAW264.7 cells stimulated with pHrodo-conjugated E. coli (20 μg/ml) showing mature endosomes. A total of 185 cells in this condition were analyzed. (B) Assessment of endosomal development at the single-cell level as determined by measuring the increase in median spot (endosome) intensity over time, which was normalized to the maximum intensities reached in that cell. (C) Single-cell endosomal maturation histograms showing the time taken to cross a threshold (as a function of threshold), which was applied to the endosomal development curves shown in (B). (D) Histograms for thresholds of 20 to 30% [that is, the first three columns of (C)] with corresponding Gaussian fits. Selected distribution is shown in red. (E) Measuring the time point of post-induction attenuation of NF-κB signaling (“off” time) of 6 randomly selected cells (from a total of 376 cells) exposed to LPS (500 ng/ml). The “off” time (determined algorithmically, not manually) for each cell is highlighted by a red circle in each trajectory. (F) Full histogram of the computed off times for all 376 measured cells treated with LPS (500 ng/ml). (G) Six randomly selected single-cell simulations of active TRIF and nuclear NF-κB trajectories in response to LPS (500 ng/ml) after introducing a cell-specific delay in endosomal maturation picked from the Gaussian distribution described in (D). (H) Predicted NF-κB signaling termination times as a function of endosomal maturation time from 200 single cells simulated as in (G). Delays of 2 hours (lower line) and 3 hours (upper) between the NF-κB signaling termination time and endosomal maturation time are indicated by red lines. All single-cell data in (A) to (F) are representative of three independent experiments.

To determine the endosomal maturation threshold that functionally inhibited TRIF signaling (Fig. 4D), we measured the timing of final NF-κB inhibition in single cells (Fig. 4E), which revealed a distribution of “off” (signal termination) times that was centered at about 6.5 to 7 hours (Fig. 4F). In our model, distributed TRIF and NF-κB signaling termination times (Fig. 4G) were closely correlated, and the effective delay between them was ~2 to 3 hours (Fig. 4H). The best-fit maturation time, then, was estimated to match a normal distribution of 4.4 ± 1.2 hours, which corresponded to an endosomal brightness threshold of 20% of the fully mature state (Fig. 4D). These results reveal an intricate relationship between constitutive and ligand-induced trafficking of TLR4 to the endosome, as well as a role for endosomal maturation in determining the transience or duration of NF-κB signaling, and they suggest that the measured cell-to-cell variability in these processes may directly contribute to the cell-to-cell variability observed in single-cell NF-κB signaling studies.

The complexity of endosomal maturation may preclude the identification of a single source of cell-to-cell variability in the process, although the maturation process is thought to be independent of the TLR signaling cascade (40). Single-cell gene expression analyses have implicated secreted interferon-β (IFN-β) in inhibiting the TLR4-dependent inflammatory response in dendritic cells (41), and we thus asked whether type I IFN receptor (IFNAR)–dependent signaling stimulated by IFN-β might accelerate endosomal maturation to inhibit TRIF signaling. However, BMDMs from IFNAR-deficient (Ifnar1−/−) mice showed identical phagosome acidification to that of wild-type cells, and spatial analysis confirmed that cells that were close to each other (that is, those cells that were exposed to similar microenvironments) were not more likely to show similar delays in endosomal maturation (fig. S9).

Extrinsic noise sources may account for the cell-to-cell variability in NF-κB dynamics

Previous work sought to explain the variability in NF-κB dynamics that occurs through intrinsic noise arising from the small number of molecules that are involved in initiating transcription of the gene encoding IκBα (42); however, in several other systems, preexisting cell-to-cell differences (extrinsic noise) were found to be the dominant determinant of variability in cellular responses and decisions (22, 43), and a sibling analysis of single-cell NF-κB dynamics supported the notion that extrinsic noise sources are also primary contributors to the variability in NF-κB dynamics (29). Similarly, we observed wide distributions in peak timing, amplitude, and duration at all concentrations of LPS tested (Fig. 5A), and we extended our modeling work to match these behaviors by focusing on extrinsic noise that affected the rates of reactions (44). Because extrinsic variation is understood to be both independent of activation and stable relative to its characteristic timescale (45), we used a distributed, deterministic model of cell-to-cell variation.

Fig. 5 The distributed, oscillatory nuclear translocation of NF-κB in response to LPS can be fitted to an extrinsic noise model.

(A) Heatmaps of measured single-cell NF-κB responses to LPS at concentrations of 5 ng/ml (left; n = 430 cells) and 5 μg/ml (right; n = 501 cells). Each row indicates the response of one cell, and the rows are ordered according to the extent of total activation. (B) The locations of reactions (indicated by number and color in the top panel) in the model topology that varied within cell population with the indicated rate distributions (bottom); these distributions represent cell-to-cell variability due to extrinsic noise: the synthesis rates of TLR4, the activation rates of TRIF and MyD88, and the endosomal maturation delay time were varied. (C and D) Five NF-κB trajectories sampled from the responses shown in (A): arrowheads in (A) indicate the trajectories of the cells shown in (C) and (D). Beside each sample trajectory, the best-fit single-cell simulations (minimum Euclidean distance) and the sampled parameter values (shown as a line indicating the specific value of the parameter relative to its parent distribution) are shown. (E) Metrics (time and amplitude of the first peak of NF-κB activity, time of the second peak, and integrated activity) were used to characterize the distributions of NF-κB responses in RAW264.7 cells (solid colors) and the best-fit simulations (gray outlines) for each of five different concentrations of LPS (500 pg/ml, 5, 50, and 500 ng/ml, and 5 μg/ml). The limits of the x and y axes were kept constant for each metric at all concentrations of LPS. (F) Fourier transformations of NF-κB dynamics for the indicated concentrations of LPS. (Left) RAW264.7 cell trajectories are shown as thin black lines, and the root mean square (RMS) of all measured responses is shown in purple. (Right) Simulations: RMS magnitude is shown in purple. Error bars show the SD from 500 simulations. All single-cell data in (A), (C), (D), and (E) are representative of three independent experiments and contain all tracked cells (n = 376 to 501 cells) in the indicated conditions.

In addition to the measured distribution of endosomal maturation delays, τmature, we identified three key receptor-proximal processes in the signal transduction cascade that might be subject to cell-to-cell differences: TLR4 synthesis, TRIF activation, and MyD88 activation (Fig. 5B). We introduced extrinsic variability into these processes by applying a zero-mean, log-normally distributed multiplier into corresponding reaction rates (Fig. 5B, bottom). These processes were considered to be representative, with log-normal distributions being the result of multiplicative combinations of multiple distributed protein species. We also considered downstream processes, such as cellular RelA abundance, but found that variability in the amount of RelA only weakly correlated with the overall dynamics of NF-κB signaling (fig. S10) (46).

A model variant that lacked endosomal maturation distribution but incorporated the three protein distributions predicted variance in the total activity of NF-κB by amplitude scaling only; knowing the amplitude of the second peak, for example, enables nearly perfect prediction of the total activity of NF-κB (r = 0.98; fig. S11A). By contrast, a model incorporating only a variable τmature predicted modulation of signal duration only; for most cells, there was no correlation between the amplitude of the second phase of NF-κB activity and its total activity (fig. S11B). The measured correlation in single-cell data between these two quantities was positive but weak, indicating combined modulation of both signal duration and amplitude (r = 0.65; fig. S11D), and the full model with varied protein and degradation delays provided a good match for this observed relationship (fig. S11C).

The full-variation model was fit to the full range of measured single-cell responses to a range of LPS concentrations (from 0.5 ng/ml to 5 μg/ml) (Fig. 5A and fig. S4). Our earlier parameterization provided us with the mean amounts of each of our distributed proteins, enabling us to fit noise levels (that is, the σ value of each log-normal distribution). To do so, we used metrics (first and second peak amplitudes, as well as total activity) to compress a large set of dynamic trajectories into univariate distributions. Model results were then optimized to fit to measured trajectories by comparing the distributions of both simulated and measured metrics, minimizing the Kolmogorov-Smirnov distance between these at all concentrations of ligand (fig. S12).

The resultant single-cell simulations qualitatively matched experimentally observed responses (Fig. 5, C and D). In response to a low concentration of LPS (5 ng/ml), almost all cells showed a first peak of NF-κB activity, with a smaller subset of cells exhibiting a second peak of NF-κB activity (Fig. 5C). At a higher concentration of LPS (5 μg/ml), almost all of the cells showed a secondary peak of NF-κB activity, with some cells continuing to exhibit oscillations in NF-κB localization for up to 8 to 9 hours (Fig. 5D). (In both cases, the values sampled from each distribution for all of the four varied parameters are shown next to the resultant simulation trajectory.) Histograms of overlaid experimental and simulation distributions of dynamical metrics showed good agreement across the full range of tested concentrations of LPS (Fig. 5E). Finally, Fourier analysis of single-cell dynamics revealed the emergence of a distinct harmonic signature, corresponding to a 2-hour periodic activity (Fig. 5F). Identical analyses performed on 500 iterations of the varied single-cell model showed a similar frequency distribution: trajectories were dominated by low-frequency information, but higher ligand concentrations gave rise to a consistent, distinct harmonic peak at roughly 0.5/hour. We note that the period encoded in these responses appears to be characteristic of the NF-κB delayed negative feedback system; similar periods have been reported in response to other stimuli (29, 47).

Variabilities in aspects of the NF-κB response are specific to either the MyD88 or TRIF pathways

The single-cell model of extrinsic noise recapitulated many aspects of the observed behavior, and could then be used to examine the individual contributions of introduced sources of noise. Simulated and measured NF-κB dynamics showed high amounts of overall activation variability, particularly at late time points (Fig. 6, A and B). We ranked both quantified and simulated dynamics by the extent of the total response and then applied the same ordering to the generated parameter values for each simulation (Fig. 6C). We showed that parameters determining TRIF activation and inactivation made the largest contributions to the response strength [correlation, r, between log(ka_TRIF) and ranking = 0.61, and correlation between τmature and ranking = 0.64], whereas ligand-TLR4 and MyD88 activation played weaker roles (r = 0.30 and 0.18, respectively).

Fig. 6 Isolating the effects of variability in specific signaling processes on the dynamics of NF-κB and IRF3.

(A) All single-cell NF-κB responses (n = 376 cells) treated with LPS (500 ng/ml) in a single representative experiment. Each row is the trajectory of one cell, and rows are ordered according to the strength of the overall response. (B) 500 single-cell model simulations of NF-κB activity, ordered according to the extent of total activity in (A), using the fully distributed parameter set (parameters were sampled from distributions in Fig. 5B). (C) The value of each of the four parameters used to generate each simulation result in rows in (B) is shown alongside that row. Each sample has been normalized to the maximum and minimum sampled values for that parameter and, except in the case of endosomal maturation, has undergone log transformation. The Spearman correlation coefficients between the sets of normalized parameter values, as well as the extent of the total NF-κB response, are shown at the top. (D) Top: Randomly generated NF-κB trajectories over time resulting from single-parameter variation models in response to LPS (500 ng/ml). Bottom: Metrics for integrated activity, amplitude and time of the first peak, amplitude of the second peak, and total duration (where [NF-κB]nucleus > 15 nM) are plotted as a function of individually varied parameters.

To examine how each noise source affected specific dynamical features of the NF-κB response, we repeated the simulations with univariate distributions (Fig. 6D, example trajectories shown on top). Variation in MyD88 activation caused modulation of the timing and amplitude of the first peak of NF-κB activity, but failed to modulate the total activity as strongly as did variations in either the activation or the inactivation of TRIF (Fig. 6D). In addition, we note that variation in MyD88 activation led to an about twofold difference in the amplitude of the first peak of NF-κB activity, which was less than the eightfold difference in the amplitude of the second peak that was caused by equivalent modulation of TRIF activation, and less than the fivefold difference in signaling duration (nuclear [NF-κB] > 0.015 μM) induced by the range of delays in endosomal maturation.

Our analysis suggests that the variability in the dynamics of NF-κB signaling in macrophages is largely dependent on variability in TRIF signaling. To better isolate the effects of this variability, we examined the dynamics of IRF3 activation, which, unlike that of NF-κB, is dependent solely on the TRIF pathway and is therefore a more direct readout of the variability in TRIF signaling. Indeed, simulating IRF3 activity in our full-variation model predicted a broadly distributed range of IRF3 responses (Fig. 7A), with a high proportion of cells showing little to no response (relative to NF-κB activity at the same concentration; compare to fig. S4; ~30% versus 5%). To test this prediction, we generated a RAW264.7 cell line expressing mVenus-IRF3, in which the nuclear translocation of IRF3 was signal-dependent. Single-cell experimental data also revealed a range of responses that were similarly distributed (Fig. 7A). We also noted unexpected oscillatory characteristics that warrant future study. The lack of a robust “on” state means that signal transduction through the TRIF pathway inherently involves information loss: no matter the magnitude of the input applied, some cells acted indistinguishably from nonactivated cells.

Fig. 7 The MyD88 and TRIF pathways encode distinct dynamics and cell-to-cell variability of NF-κB responses.

(A) Heatmap of modeled (left) and measured (right) IRF3 trajectories in 499 cells in response to LPS (50 ng/ml). Each row shows the trajectory of one cell. (B) Calculated channel capacity between outputs to a 0.5 ng/ml input (that is, no visible activation or “off”) and a varied “on” concentration of LPS. One bit indicates a perfect distinction between the “on” and “off” states in all measured single cells. Top: Simulated NF-κB dynamics in MyD88−/− and Trif−/− cells. Middle: Simulated NF-κB and IRF3 dynamics in wt cells. Bottom: NF-κB and IRF3 responses in RAW264.7 cells treated with the indicated concentrations of LPS. For IRF3 measurements, only the two highest concentrations of LPS were used. Data in (A) and (B) are representative of three independent experiments in each condition. (C) MyD88 and TRIF serve non-overlapping roles in encoding stimulus information. Given three different inputs (low, medium, and high), differential single-cell sensitivities lead to a range of encodings. Although the switch-like behavior in MyD88 activation leads to transient “all or none” responses that reliably distinguish between low and medium inputs, they fail to distinguish between medium and high. Meanwhile, TRIF encoding is highly variable, but scales upward over the three input doses. By summing the input of MyD88 and TRIF, NF-κB shows reliable levels of early activation and scalable activity afterwards.

To quantify this phenomenon, we applied information theory formalism (22) to quantify the capacities of the MyD88 and TRIF pathways to transduce information contained in stimulus inputs to the activity of transcription factor outputs. We calculated the channel capacity, a measure of the maximal fidelity of a noisy system in input-to-output information transmission (21). Specifically, we compared the output response distributions to two inputs: an “off” concentration of LPS (0.5 ng/ml) and an “on” concentration (1, 5, 50, or 500 ng/ml) (Fig. 7B). A channel capacity of 1 bit would indicate perfect “off” versus “on” encoding, enabling perfect distinction of two input conditions.

As expected, the ability to distinguish between the off and on states was improved by increasing the “on” concentration of LPS; however, our model predicted that the MyD88 pathway was less noisy than its TRIF counterpart. The NF-κB responses of both Trif−/− cells and wild-type cells showed a sharp increase in information transmission capacity when stimulated with LPS at concentrations between 1 and 5 ng/ml, or at our hypothesized “switch point” for MyD88 activation (Fig. 7B). Better information transmission was achieved for all LPS concentrations in the MyD88 pathway than when Myd88−/− cells were examined for either NF-κB or IRF3 responses, where activity was solely controlled by TRIF (Fig. 7B). Our measured data sets (1, 5, 50, and 500 ng/ml LPS for NF-κB, and 50 and 500 ng/ml LPS for IRF3) confirmed this prediction: even at high concentrations of LPS, cell-to-cell variability led to unreliable stimulus encoding of IRF3 activity (Fig. 7B, filled bars).


Here, we pursued an iterative approach of computational modeling and experimentation to develop a predictive understanding of how TLR4-responsive NF-κB dynamics in single cells are encoded by the MyD88 and TRIF pathways. We distinguished between several dynamic features and mapped these to underlying regulatory network topologies and molecular biochemical characteristics. In turn, these mechanisms determined the variability of single-cell responses in the population. Having achieved a base parameterization of the model, we added extrinsic noise to match the variability observed in the single-cell data. Specifically, we first experimentally measured one potential source of noise: endosomal maturation in single cells. We then fitted three others, related to receptor-proximal protein species, by comparing the distributions of relevant single-dimension metrics. Although this metric-match methodology could be applied to other systems, we were aided in this case by the temporal separation of our two pathways. As a result, aspects of wild-type responses mapped uniquely to the TRIF or the MyD88 pathways.

Many studies of NF-κB that address cell-to-cell variability have focused on intrinsic noise as a primary contributor (42); however, quantification in other systems (22, 28, 43, 45), and in the NF-κB system in particular (29), has suggested that the observed cell-to-cell variation originates from processes outside the TLR activation network. This extrinsic noise can be readily modeled through a distributed deterministic process, and our distributions, although constrained to four quantities, still generated a good match to the observed distributions of activity. However, we recognize that, particularly in the case of processes involving small numbers of molecules (for example, transcription or receptor activation with low amounts of ligand), intrinsic noise may also play a role in cell-to-cell variability. As we have seen here for MyD88 and TRIF, the degree to which extrinsic noise contributed to the variability of NF-κB dynamics was pathway-specific, and thus the relative contribution of intrinsic noise may also be pathway-specific. In turn, whether the appropriate modeling formalism is deterministic or stochastic will depend on the pathway under study and the question to be addressed.

We found that LPS-responsive NF-κB signaling dynamics in macrophages were oscillatory, contrary to previous observations from experiments with mouse embryonic fibroblasts (10, 13, 14). This finding was first predicted by a model of the TLR4 signaling module, which was parameterized on the basis of published quantitative data and measurements of IKK activity. It was then confirmed in single-cell experiments, which also revealed substantial cell-to-cell variability that explained why such oscillatory dynamics were not evident in previous measurements at the population level. Thus, oscillatory dynamics may be a more profound and conserved feature of NF-κB signaling than was previously thought, and they may be independent of the stimulus and solely a property of IκB feedback, which was previously derived theoretically (47).

Structural characterization of MyD88 (37) revealed that it oligomerizes into a large signaling molecular complex or signalosome. This oligomerization was hypothesized to generate positive cooperativity (48, 49), and similar signalosome-based cooperativity is thought to generate bistability in the decision-making process for apoptosis (50). However, although signalosome-forming adaptor proteins provide high specificity by selectively recruiting substrates (51), their effects on the dynamics of signal transduction remain unclear (37, 52). Our study suggests that MyD88 signalosome–mediated cooperativity, modeled by and fitted to a Hill equation, is sufficient to explain a concentration threshold in both IKK activity and NF-κB translocation. This cooperativity ensured not only the reliability of MyD88 signaling but also its transience (Fig. 2G, together with the translocation module, and Fig. 3) while reducing its scalability in response to different ligand concentrations (Fig. 2H). Several other signal transmission systems upstream of NF-κB exhibit similar thresholding behavior (which is mediated by either cooperativity or positive feedback), including in the TNF receptor 1–IKK signal transduction axis (53) as well as in B cell receptor (BCR) signaling (54), in which the scaffold protein CARMA1 mediates IKK-dependent positive feedback. Note that the positive feedback motif in BCR signaling does not ensure the transience, but rather the prolongation, of NF-κB signaling dynamics.

Receptor endocytosis provides temporal separation between NF-κB activation by the MyD88 and TRIF pathways. Although this delay is small relative to the oscillatory period caused by IκB-dependent feedback, it is sufficiently long to ensure that the timing of the first peak of NF-κB activity is largely driven by MyD88 signaling, which thereby determines the oscillatory phase of the cell. Because LPS, CD14, MD2, and TLR4 form stable complexes (55), and because LPS itself is stable, receptor trafficking serves to determine the duration of MyD88 signaling, and endosomal maturation serves to determine the duration of TRIF signaling. Our work suggests that the process of deactivation of TLR4 through endosome maturation plays a role in determining the total duration of signaling. Single-cell measurements of this process showed high variability, adding another source of noise in encoding late-phase NF-κB activity. Inflammatory gene expression programs induced by TLR4 signaling reach their maximal induction at different times (56, 57) and may differentially depend on NF-κB signaling dynamics. Our study suggests that genes dependent on the first peak of NF-κB activity for their expression may show less cell-to-cell variability than those dependent on late-phase NF-κB activity. Consistent with this suggestion, a single-cell study of gene expression noted subsets of genes that were unimodally, universally expressed early in dendritic cells that were exposed to LPS, which was followed by a multimodal, highly variable distribution of inflammatory gene expression at later times (41).

Innate immune responses must ensure both the high sensitivity and tight control that provide appropriate responses that minimize the potentially destructive consequences to surrounding tissues (58). The integration of two pathways with distinct signaling topologies and system characteristics may enable the TLR4 network to accomplish these competing tasks. Our results suggest that in the face of high extrinsic noise associated with macrophage and microenvironment heterogeneity, the MyD88 pathway is still able to provide a reliable first signaling response because of signalosome topology (Fig. 7C, left); however, this is limited in duration and not scalable with ligand concentration, thus minimizing the risk of inflammatory damage. In contrast, the TRIF pathway exploits cell-to-cell variability to provide a prolonged response only in a fraction of cells, and this response is scalable with ligand concentration (Fig. 7C, right). We speculate that by limiting the number of cells that can produce cytokines at amounts sufficient to potentially act systemically, the network behavior also reduces the risk of inflammatory shock while providing for scalability in systemic immune activation. Work delineating the control of expression and functions of TNF is consistent with this view. Whereas local autocrine functions of TNF only require MyD88-mediated signaling events, TRIF-mediated signaling mechanisms are required for the full production of TNF and its (nonlocal) paracrine functions (16).

By combining the two pathway responses (Fig. 7C, bottom middle), every cell in a population is able to mount a minimal response to LPS to provide a digital and reliable local response, and this response is able to scale (in subsets of cells) with the extent of the threat, tuning systemic immune functions appropriately. Our study thus reveals how distinct topologies of innate immune signaling pathways either mitigate or exploit molecular network noise to provide appropriate responses at the cell, tissue, or organism level that are appropriate to the pathogen threat while carefully controlling the risk of endotoxic shock.


Mathematical model

Three modules constitute the model: the TLR4 module, the IRF3 module, and the NF-κB module (Fig. 1A). The details of the model, which is modeled as chemical reaction networks (see also table S1), are shown in fig. S2. The IKK and NF-κB modules (parameters labeled by green numbers in fig. S2) are adapted from our previous studies (33, 34). All of the components of the TLR4 and IRF3 modules were newly formalized by ODEs based on mass action, except for the case of MyD88 activation, for which the Hill equation was used to represent Myddosome formation. We parameterized the ligand-receptor binding and shuttling parts of the model with data from the published literature (parameters labeled 1 to 15 in fig. S2 and Supplementary text). The model is numerically implemented in MATLAB R2014a (MathWorks), and the remaining parameters are fitted with lsqnonlin, a constrained minimization algorithm from the MATLAB Optimization Toolbox (MathWorks), using the quantified BMDM IKK and IRF3 data. The final parameter values, together with the chemical reactions, are shown in table S1.

Simulating single-cell trajectories in RAW264.7 cells

To parameterize the model to single-cell data, we first accounted for sensitivity differences between BMDMs and RAW264.7 cells by multiplying the abundances of MyD88 and TRIF in BMDMs by 0.33 and 0.34, respectively. Next, we added preexisting cell-to-cell variability (extrinsic noise, μ = 0 and fitted σ for each log-normal distribution) at discrete points in the model: upstream (TLR4 synthesis), the MyD88 branch (activation), and the TRIF branch (activation). Finally, we explicitly modeled variable termination of TLR4 signaling as a result of endosomal maturation, which was parameterized by single-cell measurements that revealed timing differences in the initiation of maturation (Fig. 4). These distributions were used to guide a normally distributed delay parameter (τmature), after which endosomal TLR4-LPS was rapidly degraded (with a rate constant of 2 min–1). The simulated single-cell trajectories in Fig. 5, C and D, were selected based on minimizing the Gaussian distance between the simulation and experimental trajectories.

Cell culture and biochemical assays

Wild-type, Myd88−/−, and Trif−/− C57BL/6 mice were housed at the University of California, San Diego, in accordance with protocols authorized by the Institutional Animal Care and Use Committee. BMDMs were generated by culturing 6 × 106 bone marrow cells from mouse femurs in suspension in L929-conditioned Dulbecco’s modified Eagle’s medium (DMEM) on 15-cm plates 7 days before replating them for experiments to be performed on day 8. LPS was sourced from Sigma (B5:055). Antibodies against RelA/p65 (sc-372) were from Santa Cruz Biotechnology. EMSA and kinase assays were performed as previously described (59, 60). RAW264.7 cells (ATCC TIB-71) were sequentially transduced with two lentiviral vectors encoding EYFP-RelA (under the control of the 3.5-kb sequence upstream of native RelA) and H2B-mCherry. Double-stable lines were made by successive selection and then were further purified by fluorescence-activated cell sorting. Only cells from passages 16 to 20 were used for imaging experiments. RAW264.7 cells were maintained in DMEM (CellGro 10-013) supplemented with 10% fetal bovine serum, 20 mM Hepes, and 1× penicillin/streptomycin. Twenty hours before starting the experiments, the cells were replated in eight-well μ-slides (ibidi) at a density of 50,000 cells/cm2. Two hours before the experiment began, one-third of the total volume of medium was drawn off and mixed with the appropriate stimulus, which was then injected into the chamber precisely at the start of the experiment.


Cells, plated on eight-well μ-slides (ibidi), were incubated at 5% CO2 and 37°C and imaged every 5 min with a Zeiss AxioObserver fitted with a 40× oil immersion objective, LED (light-emitting diode) fluorescence excitation, and CoolSnap HQ2 camera. We collected differential interference contrast, mCherry, and YFP images over 18.5 hours at 12 stage positions per experimental condition, and then exported raw images to MATLAB for automated analysis, which was performed as previously described (22).


Supplementary Text

Fig. S1. Fitting the model to biochemical measurements of BMDMs.

Fig. S2. Full schematic diagram of the TLR4–NF-κB model reaction network.

Fig. S3. Feed-forward loops from secreted TNF and IL-1 do not affect the dynamics of NF-κB signaling in response to LPS.

Fig. S4. Full range of the measured, single-cell NF-κB responses from live-cell imaging experiments.

Fig. S5. The model of Myddosome formation can be approached by Hill kinetics.

Fig. S6. Assay of the LPS concentration–dependent stimulation of the kinase activity of IKK.

Fig. S7. Robustness of the primary shuttling mode to parameter perturbation.

Fig. S8. Fitting and modeling of the endosomal maturation–induced degradation of TLR4 complexes.

Fig. S9. Endosomal maturation does not appear to be controlled by paracrine IFN signaling.

Fig. S10. RelA abundance is unrelated to the dynamics of NF-κB signaling.

Fig. S11. Predictions of the dynamics of NF-κB signaling based on varied degradation rates and protein abundances.

Fig. S12. Fitting extrinsic variability to aspects of NF-κB dynamics.

Table S1. Model reactions and parameter values.

Computational Model Files

References (6174)


Acknowledgments: We thank S. Mitchell and J. Vargas for proofreading and guidance. Funding: The work was supported by NIH grants R01 GM071573 and R01 CA141722. Author contributions: Z.C. developed the mathematical model with contributions from B.T.; B.T. generated all of the single-cell data; D.R.O. contributed biochemical measurements; and Z.C., B.T., and A.H. wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Model code is available at; images and other materials are available from A.H.
View Abstract

Stay Connected to Science Signaling

Navigate This Article