Biochemical feedback and its application to immune cells II: dynamics and critical slowing down

TL;DR

This paper explores how biochemical feedback in immune cells exhibits critical slowing down near bifurcation points, using a mapping to the Ising model to quantify response times and hysteresis effects during perturbations.

Contribution

It introduces a novel application of the Ising universality class to model critical slowing down in biochemical feedback systems, linking theory to experimental immune cell responses.

Findings

Critical slowing down observed near bifurcation points in immune cell models.

Quenches closer to the critical point lead to slower cellular responses.

The Kibble-Zurek mechanism predicts hysteresis scaling in cellular responses.

Abstract

Near a bifurcation point, the response time of a system is expected to diverge due to the phenomenon of critical slowing down. We investigate critical slowing down in well-mixed stochastic models of biochemical feedback by exploiting a mapping to the mean-field Ising universality class. This mapping allows us to quantify critical slowing down in experiments where we measure the response of T cells to drugs. Specifically, the addition of a drug is equivalent to a sudden quench in parameter space, and we find that quenches that take the cell closer to its critical point result in slower responses. We further demonstrate that our class of biochemical feedback models exhibits the Kibble-Zurek collapse for continuously driven systems, which predicts the scaling of hysteresis in cellular responses to more gradual perturbations. We discuss the implications of our results in terms of the tradeoff between a precise and a fast response.

Figures

Equations

$$j_{r}X+Y_r^{+}\rightleftharpoons (j_r+1)X+Y_r^{-},$$

$$p_n=\frac{p_0}{n!}\prod _{j=1}^n{f}_j,$$

$$f_n=\frac{n{p}_n}{p_{n-1}},$$

$$m\equiv \frac{n_{\ast }-n_c}{n_c},h\equiv \frac{2(f_{n_c}-n_c)}{-f_{n_c}^{\prime \prime \prime }{n}_c^3},\theta \equiv \frac{2(1-f_{n_c}^{\prime })}{-f_{n_c}^{\prime \prime \prime }{n}_c^2},$$

$$A\underset{k_1^{-}}{\overset{k_1^{+}}{\rightleftharpoons }}X,2X+B\underset{k_2^{-}}{\overset{k_2^{+}}{\rightleftharpoons }}3X.$$

$$f_n=\frac{a{K}^2+s(n-1)(n-2)}{(n-1)(n-2)+K^2},$$

$${\tau }_c{|}_{h=0}$$

$${\tau }_c{|}_{\theta =0}$$

$${\tau }_c=\frac{1}{\kappa (0)}{\int }_0^{\infty }dt\kappa (t),$$

$${\tau }_r=\frac{1}{\Delta \bar{n}(0)}{\int }_0^{t_{\max }}dt\Delta \bar{n}(t),$$

$$d_c={\left[({\theta }^{\beta \delta }{)}^2+h^2\right]}^{1/2}.$$

$${\tau }_r=\frac{1}{\Delta \bar{I}(0)}{\int }_0^{t_{\max }}dt\Delta \bar{I}(t),$$

$$\theta$$

$$h$$

$$m$$

$$n_c$$

$$\dot{\vec{p}}=\mathbf{L}\vec{p}.$$

$$\mathbf{L}{\vec{v}}_j$$

$${\vec{u}}_j\mathbf{L}$$

$$p_n(t)=\sum _{j{n}^{\prime }}{u}_{j{n}^{\prime }}{p}_{n^{\prime }}(0)e^{{\lambda }_{j}t}{v}_{nj}$$

$$\kappa (t)=-{\bar{n}}^2+\sum _{mn}{p}_{mn}mn=-{\bar{n}}^2+\sum _{mn}{p}_{n|m}{p}_{m}mn,$$

$$\kappa (t)$$

$$=\sum _{mn}mn{v}_{m0}\sum _{j=1}^{\infty }{u}_{jm}{e}^{{\lambda }_{j}t}{v}_{nj},$$

$${\tau }_c=\frac{1}{{\sigma }^2}\sum _{mn}mn{v}_{m0}\sum _{j=1}^{\infty }{u}_{jm}\left(\frac{1}{-{\lambda }_j}\right)v_{nj}.$$

$${\tau }_c={\sigma }^{-2}\vec{n}\mathbf{V}\mathbf{F}\mathbf{U}\vec{w},$$

$${\tau }_c=\frac{{\tau }_b{\sigma }_b^2}{2{\sigma }^2},$$

Full text

I Introduction

Critical slowing down is the phenomenon in which the relaxation time of a dynamical system diverges at a bifurcation point Strogatz (2018). Biological systems are inherently dynamic, and therefore one generally expects critical slowing down to accompany transitions between their dynamic regimes. Indeed, signatures of critical slowing down, including increased autocorrelation time and increased fluctuations, have been shown to precede an extinction transition in many biological populations Scheffer et al. (2009, 2012), including bacteria Veraart et al. (2012), yeast Dai et al. (2012), and entire ecosystems Wang et al. (2012). Similar signatures are also found in other biological time series, including dynamics of protein activity Sha et al. (2003) and neural spike dynamics Meisel et al. (2015).

Canonically, critical slowing down depends on scaling exponents that define divergences along particular parameter directions in the vicinity of a critical point Hohenberg and Halperin (1977). Therefore, connecting the theory of critical slowing down to biological data requires identification of thermodynamic state variables, their scaling exponents, and a principled definition of distance from the critical point. However, in most biological systems it is not obvious how to define the thermodynamic state variables, let alone scaling exponents and distance from criticality. In a previous study Erez et al. (arXiv:1703.04194) we showed how near its bifurcation point, a class of biochemical systems can be mapped to the mean-field Ising model, thus defining the state variables and their associated scaling exponents. This provides a starting point for the investigation of critical slowing down in such systems, as well as how to apply such a theory to experimental data.

Additionally, most studies of critical slowing down in biological systems investigate the response to a sudden experimental perturbation (a “quench”), such as a dilution or the addition of a nutrient or drug. This leaves unexplored the response to gradual environmental changes, a common natural scenario. When a gradual change drives a system near its critical point, critical slowing down delays the system’s response such that no matter how gradual the change, the response lags behind the driving. In physical systems this effect is known as the Kibble-Zurek mechanism Kibble (1976); Zurek (1985), which predicts these nonequilibrium lagging dynamics in terms of the exponents of the critical point. It remains unclear whether and how the Kibble-Zurek mechanism applies to biological systems.

Here we investigate critical slowing down for well-mixed biochemical networks with positive feedback, and we use our theory to interpret the response of immune cells to an inhibitory drug. Using our previously derived mapping Erez et al. (arXiv:1703.04194), we show theoretically that critical slowing down in our class of models proceeds according to the static and dynamic exponents of the mean-field Ising universality class. The mapping identifies an effective temperature and magnetic field in terms of the biochemical parameters, which defines a distance from the critical point that can be extracted from experimental fluorescence data. We find that drug-induced quenches that take an immune cell closer to its critical point result in longer response times, in qualitative agreement with our theory. We then show theoretically that our system, when driven across its bifurcation point, falls out of steady state in the manner predicted by the Kibble-Zurek mechanism, thereby extending Kibble-Zurek theory to a biologically relevant nonequilibrium setting. Our work elucidates the effects of critical slowing down in biological systems with feedback, and provides insights for interpreting cell responses near a dynamical transition point.

II Results

We consider a well-mixed reaction network in a cell where $X$ is the molecular species of interest, and the other species $A$, $B$, $C$, etc. form a chemical bath for $X$ [Fig. 1(a)]. Whereas previously we considered only the steady state distribution of $X$ Erez et al. (arXiv:1703.04194), here we focus on dynamics in and out of steady state. Specifically, as shown in Fig. 1(b), we consider (i) steady state, where the bath is constant in time; (ii) a quench, where the bath changes its parameters suddenly; and (iii) driving, where the bath changes its parameters slowly and continuously. In each case we are interested in a corresponding timescale: (i) the autocorrelation time $\tau_{c}$ of $X$, (ii) the response time $\tau_{r}$ of $X$, and (iii) the driving time $\tau_{d}$ of the bath.

First we review the key features of our stochastic framework for biochemical feedback and its mapping to the mean-field Ising model Erez et al. (arXiv:1703.04194). We consider an arbitrary number of reactions $r$ in which $X$ is produced from bath species $Y_{r}^{\pm}$ and/or $X$ itself (feedback),

[TABLE]

where $j_{r}$ are stoichiometric integers. The probability of observing $n$ molecules of species $X$ in steady state according to Eq. 1 is

[TABLE]

where $p_{0}^{-1}=\sum_{n=0}^{\infty}(1/n!)\prod_{j=1}^{n}f_{j}$ is set by normalization, and $f_{n}$ is a nonlinear feedback function governed by the reaction network. The inverse of Eq. 2,

[TABLE]

allows calculation of the feedback function from the distribution. The function $f_{n}$ determines an effective order parameter, reduced temperature, and magnetic field,

[TABLE]

respectively, where $n_{c}$ is defined by $f^{\prime\prime}_{n_{c}}=0$, and $n_{*}$ are the maxima of $p_{n}$. Qualitatively, $n_{c}$ sets the typical molecule number, $\theta$ drives the system to a unimodal ($\theta>0$) or bimodal ($\theta<0$) state, and $h$ biases the system to high ($h>0$) or low ($h<0$) molecule numbers. The critical point occurs at $\theta=h=0$. The state variables $m$, $\theta$, and $h$ scale according to the exponents $\alpha=0$, $\beta=1/2$, $\gamma=1$, and $\delta=3$ of the mean-field Ising universality class. Detailed analysis of this mapping in steady state is found in our previous work Erez et al. (arXiv:1703.04194).

Near the critical point, all specific realizations of a class of systems scale in the same way, and therefore it suffices to consider a particular realization of Eq. 1 from here on. We choose Schlögl’s second model Erez et al. (arXiv:1703.04194), a simple and well-studied case Schlögl (1972); Dewel et al. (1977); Nicolis and Malek-Mansour (1980); Brachet and Tirapegui (1981); Grassberger (1982); Prakash and Nicolis (1997); Liu et al. (2007); Vellela and Qian (2009) in which $X$ is either produced spontaneously from bath species $A$, or in a trimolecular reaction from two existing $X$ molecules and bath species $B$,

[TABLE]

In this case the feedback function is

[TABLE]

where we have introduced the dimensionless quantities $a\equiv k_{1}^{+}n_{A}/k_{1}^{-}$, $s\equiv k_{2}^{+}n_{B}/k_{2}^{-}$, and $K^{2}\equiv k_{1}^{-}/k_{2}^{-}$ in terms of the reaction rates and the numbers of $A$ and $B$ molecules. Given Eqs. 4 and 6, the effective thermodynamic variables $n_{c}$, $\theta$, and $h$ can be written in terms of $a$, $s$, and $K$ or vice versa Erez et al. (arXiv:1703.04194), with $1/k_{1}^{-}$ setting the units of time.

II.1 Critical slowing down in steady state

In steady state, critical slowing down causes correlations to become long-lived near a dynamical transition point. Qualitatively, the fixed point is transitioning from stable to unstable, and therefore the basin of attraction is becoming increasingly wide. As a result, a dynamic trajectory takes increasingly long excursions from the mean, making it heavily autocorrelated. The autocorrelation time $\tau_{c}$ diverges at the critical point according to Pathria and Beale (2011)

[TABLE]

where we expect $\nu z=1$ for mean-field dynamics Hohenberg and Halperin (1977); Kopietz et al. (2010). Here the autocorrelation time $\tau_{c}$ is defined as

[TABLE]

where $\kappa(t)=\langle n(0)n(t)\rangle-\bar{n}^{2}$ is the steady-state autocorrelation function, $\kappa(0)=\sigma^{2}$ is the variance, and we have taken the start time to be $t=0$ without loss of generality because the system is in steady state.

To confirm the value of $\nu z$, we plot $\tau_{c}$ vs. $h$ at $\theta=0$ (Eq. 8). We calculate $\tau_{c}$ either directly from the master equation or from stochastic simulations Gillespie (1977) using the method of batch means Thompson (2010) (see Appendix A). The results are shown in Fig. 2. We see in Fig. 2(a) that $\tau_{c}$ indeed diverges with $h$, and that the location of the divergence approaches the expected value $h=0$ as the molecule number $n_{c}$ increases. We also see that the height of the peak increases with $n_{c}$ due to the rounding of the divergence Stephens et al. (2013). The inset of Fig. 2(b) plots this dependence: we see that $\tau_{c}$ at the critical point $\theta=h=0$ scales like $n_{c}^{1/2}$ for large $n_{c}$ (the application of this dependence to dynamic driving will be discussed in Section II.3). Finally, we see in the main panel of Fig. 2(b) that when $n_{c}$ is sufficiently large, $\tau_{c}$ falls off with $|h|$ with the expected scaling exponent of $\nu z/\beta\delta=2/3$. Taken together, these results confirm that the divergence of the autocorrelation time in the Schlögl model obeys the static exponents of the mean-field Ising universality class ($\beta\delta=3/2$) and the dynamic expectation for mean-field systems ($\nu z=1$).

II.2 Quench response and application to immune cells

When subjected to a sudden environmental change (a quench), the system will take some finite amount of time to respond [Fig. 1(b), middle]. We expect that if a quench takes the system closer to its critical point, the response time should be longer due to critical slowing down. To make this expectation quantitative, we define the response time $\tau_{r}$ in terms of the dynamics of the mean molecule number $\bar{n}$ as

[TABLE]

where the quench occurs at $t=0$, we define $\Delta\bar{n}(t)=\bar{n}(t)-\bar{n}(t_{\max})$, and we ensure that $t_{\max}\gg\tau_{r}$. We compute $\bar{n}(t)$ from the time-dependent distribution $p_{n}(t)$ using the stochastic simulations. Examples of $p_{n}(t)$ for a small and a large quench are shown in Fig. 3(a).

We define the distance from the critical point in terms of the state variables $\theta$ and $h$. Specifically, $\tau_{c}$ scales identically with $\theta^{\beta\delta}$ as it does with $h$ (Eqs. 7 and 8), which defines the Euclidean distance $d_{c}$ from the critical point as

[TABLE]

This measure will be important when comparing with the experiments because, as opposed to in most condensed matter experiments, it is difficult in the biological experiments we describe to manipulate only one parameter ($\theta$ or $h$) independently of the other.

To test whether the response time increases with proximity to the critical point, we must define initial values $\theta_{0}$ and $h_{0}$ for the environment before the quench, and a series of values $\theta$ and $h$ for the environment after the quench that are varying distances from the critical point $\theta=h=0$. There are many such choices for these values, but anticipating the experimental results that we will describe shortly, we choose the initial point (black square) and final points (colored circles) shown in Fig. 3(b). Dotted curves of equal $d_{c}$ are also shown, which make clear that larger quenches (yellow circles) take the system farther from the critical point than smaller quenches (blue circles). The dependence of $\tau_{r}$ on $d_{c}$ is shown in Fig. 3(c), and we see that indeed $\tau_{r}$ decreases as $d_{c}$ increases, or equivalently the response time increases with proximity to the critical point.

We now compare our theory with data from immune cells. We focus on the abundance in T cells of doubly phosphorylated ERK (ppERK), a protein that initiates cell proliferation and is implicated in the self/non-self decision between mounting an immune response or not Vogel et al. (2016); Altan-Bonnet and Germain (2005). Specifically, we use flow cytometry to measure the ppERK distribution at various times after the addition of a drug that inhibits SRC, a key enzyme in the cascade that leads to ERK phosphorlyation (see Appendix B for experimental methods). When the dose of the drug is small, the distribution hardly changes [Fig. 3(d), top]; whereas when the dose is large, the distribution changes significantly [Fig. 3(d), bottom]. The responses to all doses are shown in Appendix B.

After the addition of the drug, the cells reach a new steady-state ppERK distribution [green curves in Fig. 3(d)]. The distribution corresponds to an effective feedback function via Eq. 3, from which the effective temperature $\theta$ and field $h$ can be calculated via Eq. 4 Erez et al. (arXiv:1703.04194). The values of $\theta$ and $h$ calculated from the experimental distributions are shown in Fig. 3(e). We see that larger doses take the cells farther from their initial distribution (black square), as expected. We also see that larger doses take the system farther from the critical point $\theta=h=0$. The general shape of the $\theta$ and $h$ values motivated our choice of theoretical values in Fig. 3(b).

We define the response time to the drug as in Eq. 10, here in terms of the mean fluorescence intensity of ppERK,

[TABLE]

where $\Delta\bar{I}(t)=\bar{I}(t)-\bar{I}(t_{\max})$ and $t_{\max}=30$ min. We calculate the distance from criticality using Eq. 11 as before, here using the experimental values of $\theta$ and $h$. We see in Fig. 3(f) that the response time $\tau_{r}$ decreases with the distance from criticality $d_{c}$, consistent with the prediction from the theory [Fig. 3(c)]. This suggests that critical slowing down occurs in the response of the T cells to the drug.

Although in Fig. 3(f) the response time $\tau_{r}$ comes directly from the experimental data, the distance from criticality $d_{c}$ is calculated from the experimental data using expressions from the theory (Eqs. 3 and 4). This makes the results in Figs. 3(c) and 3(f) not entirely independent. To confirm that the agreement between Figs. 3(c) and 3(f) is not a result of an implicit co-dependence, we seek a measure that is related to distance from criticality but that is not dependent on the theory. We choose the entropy of the distribution $S=-\sum_{n}p_{n}\log p_{n}$ because near criticality, the distribution is broad and flat, and therefore we expect the entropy to be large; whereas far from criticality, the distribution has either one or two narrow peaks, and therefore we expect the entropy to be small Erez et al. (arXiv:1703.04194). Indeed, we see in Fig. 3(g) that in the theory, the response time $\tau_{r}$ increases with the entropy $S$, consistent with the fact that it decreases with the distance from criticality [Fig. 3(c)]. The same is evident in the experiments: we see in Fig. 3(h) that low drug doses correspond to long response times and high entropies, whereas high drug doses correspond to short response times and low entropies, resulting in an increase of response time $\tau_{r}$ with entropy $S$. Calculating the entropy in Fig. 3(h) requires a conversion between intensity $I$ and molecule number $n$, and we have checked that the results in Fig. 3(h) are qualitatively unchanged for different choices of this conversion factor over several orders of magnitude. The agreement between Figs. 3(g) and 3(h) offers further evidence that the T cells experience critical slowing down, with the data analysis completely independent from our theory.

II.3 Dynamic driving and Kibble-Zurek collapse

While some environmental changes are sudden, many changes in a biological context are gradual [Fig. 1(b), right]. When a gradual change drives a system through its critical point, critical slowing down delays the system’s response such that no matter how gradual the change, the response lags behind the driving. Although in a biological setting the driving protocol could take many forms, terms beyond the leading-order linear term do not change the critical dynamics Chandran et al. (2012). This is a major theoretical advantage because it allows us to specialize to linear driving without loss of biological realism. Specifically, we focus on linear driving across the critical point with driving time $\tau_{d}$, setting either $\theta(t)=\theta_{i}-(\theta_{f}-\theta_{i})t/\tau_{d}$ and $h=0$, or $h(t)=h_{i}-(h_{f}-h_{i})t/\tau_{d}$ and $\theta=0$, where $i$ and $f$ denote the initial and final parameter values, respectively.

In a traditional equilibrium setting, the dynamics of lagging trajectories are described in terms of the critical exponents by the Kibble-Zurek mechanism Kibble (1976); Zurek (1985). The idea of the Kibble-Zurek mechanism is that far from the critical point, the change in the system’s correlation time due to the driving, over a correlation time, is small compared to the correlation time itself, $(d\tau_{c}/dt)\tau_{c}\ll\tau_{c}$, and therefore the system responds adiabatically. However, as the system is driven closer to the critical point, these two quantities are on the same order, or $d\tau_{c}/dt\sim 1$, and the system begins to lag. Applying this condition to Eqs. 7 and 8, and using the above expressions for $\theta(t)$ and $h(t)$, one obtains

[TABLE]

respectively. Because $m\sim(-\theta)^{\beta}$ or $m\sim h^{1/\delta}$ near criticality in the mean-field Ising class, we have

[TABLE]

respectively. Therefore, if the system is driven at different timescales $\tau_{d}$, the Kibble-Zurek mechanism predicts that plots of the rescaled variables $m\tau_{d}^{\beta/(\nu z+1)}$ vs. $\theta\tau_{d}^{1/(\nu z+1)}$ or $m\tau_{d}^{\beta/(\nu z+\beta\delta)}$ vs. $h\tau_{d}^{\beta\delta/(\nu z+\beta\delta)}$ will collapse onto universal curves.

When testing these predictions using simulations of a spatially extended physical system, the finite size of the system causes a truncation of the autocorrelation time. This truncation is usually accounted for using a finite-size correction Chandran et al. (2012). In our system, a similar truncation of the autocorrelation time is caused by the finite number of molecules. Specifically, the inset of Fig. 2(b) shows that at criticality we have $\tau_{c}\sim n_{c}^{1/2}$ for large $n_{c}$, where $n_{c}$ sets the typical number of molecules in the system. Therefore, we interpret $n_{c}$ as a “system size,” and we correct for finite-size effects in the following way. Combining the relation $\tau_{c}\sim n_{c}^{1/2}$ with Eqs. 7 and 8, and Eqs. 13 and 14, we obtain

[TABLE]

for the driving of $\theta$ or $h$, respectively. We choose $n_{c}$ arbitrarily for a particular driving time $\tau_{d}$, and when we choose a new $\tau_{d}$, we scale $n_{c}$ appropriately according to Eqs. 17 and 18.

This procedure allows us to test the predictions of the Kibble-Zurek mechanism using simulations of the Schlögl model. The results are shown in Fig. 4. We see in Fig. 4(a) that as $\theta$ is driven from a positive to a negative value, the bifurcation response is lagging, occurring at a value less than the critical value $\theta=0$ (supercooling). Conversely, when $\theta$ is driven from a negative to a positive value, the convergence occurs at a value greater than $\theta=0$ (superheating). In both directions, the lag is larger when the driving is faster, corresponding to smaller values of $\tau_{d}$ (from yellow to dark brown). We see in Fig. 4(b) that similar effects occur for the driving of $h$. Yet, we see in Figs. 4(c) and (d) that the rescaled variables collapse onto single, direction-dependent curves within large regions near criticality. Note that the direction dependence (i.e., hysteresis) is preserved as part of these universal curves, but the lags vanish in the collapse. This result demonstrates that our nonequilibrium birth-death model exhibits the Kibble-Zurek collapse predicted for critical systems. Together with our previous findings, this result suggests that such a collapse should emerge in biological experiments where environmental parameters (e.g., drug dose) are dynamically controlled in a gradual manner. More broadly, by phenomenologically collapsing such experimental curves, it should be possible to deduce the critical exponents of such biological systems without fine-tuning them to criticality, but instead by gradual parameter sweeps.

III Discussion

We have investigated critical slowing down in a minimal stochastic model of biochemical feedback. By exploiting a mapping to Ising-like thermodynamic variables, we have made quantitative predictions for the response of a system with feedback to both sudden and gradual environmental changes. In response to a sudden change (a quench), we have shown that the system will respond more slowly if the quench takes it closer to its critical point, in qualitative agreement with multiple-time-point flow cytometry experiments in immune cells. In response to more gradual driving, we have shown that the lagging dynamics of the system proceed according to the Kibble-Zurek mechanism for driven critical phenomena. Together, our results elucidate the consequences of critical slowing down for biochemical systems with feedback, and demonstrate those consequences on an example system from immunology.

For the immune cells, critical slowing down may present a tradeoff in terms of the speed vs. the precision of an immune response. Specifically, ppERK is implicated in the decision of whether or not to mount the immune response Vogel et al. (2016); Altan-Bonnet and Germain (2005), suggesting that ppERK dynamics near the bifurcation point are of key biological importance. Yet, the bifurcation point is the point where critical slowing down is most pronounced. In fact, the inset of Fig. 2(b) demonstrates that the system slows down as the number of molecules in the system increases. On the other hand, large molecule number is known to decrease intrinsic noise and thereby increase the precision of a response Elowitz et al. (2002). This suggests that cells may face a tradeoff in terms of speed vs. precision when responding to changes that occur near criticality, as suggested for other biological systems Skoge et al. (2011); Mora and Bialek (2011).

Our work extends the Kibble-Zurek mechanism to a nonequilibrium biological context. Traditionally, the mechanism has been applied to physical systems from cosmology Kibble (1976) and from hard Zurek (1985); Del Campo and Zurek (2014) or soft Deutschländer et al. (2015) condensed matter. Here, we extend the mechanism to the context of biochemical networks with feedback, where the system already exists in a nonequilibrium steady state, and the external protocol takes the system further out of equilibrium into a driven state. It will be interesting to see to what other nonequilibrium contexts the Kibble-Zurek mechanism can be successfully applied Deffner (2017).

The theory we present here assumes only intrinsic birth-death reactions and neglects more complex mechanisms such as bursting Friedman et al. (2006); Mugler et al. (2009), parameter fluctuations Shahrezaei et al. (2008); Horsthemke and Lefever (1984), or cell-to-cell variability Cotari et al. (2013); Erez et al. (2018) that may play an important role in the immune cells. Nonetheless, similar models that also focus only on intrinsic noise have successfully described ppERK in T cells in the past Das et al. (2009); Prill et al. (2015). Moreover, we expect that intrinsic fluctuations should play their largest role near the bifurcation point. Finally, we expect that near the bifurcation point, the essential behavior of the system should be captured by any model that falls within the appropriate universality class.

In this and previous work Erez et al. (arXiv:1703.04194) we have explored the dynamic and static scaling properties of single cells subject to biochemical feedback. Natural extensions include generalizing the theory to cell populations or other systems that are not well-mixed such as intracellular compartments. This would allow one to investigate the spatial consequences of proximity to a bifurcation point, such as long-range correlations in molecule numbers and the associated implications for sensing, information transmission, patterning, or other biological functions.

Acknowledgments

We thank Anushya Chandran for helpful communications. This work was supported by Simons Foundation grant 376198 (T.A.B. and A.M.), Human Frontier Science Program grant LT000123/2014 (Amir Erez), National Science Foundation Research Experiences for Undergraduates grant PHY-1460899 (C.P.), National Institutes of Health (NIH) grants R01 GM082938 (A.E.) and R01 AI083408 (A.E., R.V., and G.A.-B) and the NIH National Cancer Institute Intramural Research programs of the Center for Cancer Research (A.E. and G.A.-B.).

Appendix A Autocorrelation time

We calculate the autocorrelation time $\tau_{c}$ (Eq. 9) for the Schlögl model in steady state using one of two methods, the first more efficient for small molecule numbers, and the second more efficient for large molecule numbers. The first method is to calculate $\tau_{c}$ numerically from the master equation for $p_{n}$ by eigenfunction expansion. The master equation follows from the reactions in Eq. 5 Erez et al. (arXiv:1703.04194) and can be written in vector notation as

[TABLE]

where ${\bf L}$ is a tridiagonal matrix containing the birth and death propensities for $X$. The eigenvectors of L satisfy

[TABLE]

where the eigenvalues obey $\lambda_{j}\leq 0$ with only $\lambda_{0}$ vanishing for the steady state, and $\vec{v}_{j}^{T}\neq\vec{u}_{j}$ because L is not Hermitian Walczak et al. (2009). Because Eq. 19 is linear in $\vec{p}$, the solution is

[TABLE]

for initial condition $p_{n}(0)$. Calling $n(0)\equiv m$ and $n(t)\equiv n$, we write the autocorrelation function (see Eq. 9) as

[TABLE]

where $p_{m}=v_{m0}$ is the steady-state distribution, and $p_{n|m}$ is the dynamic solution at time $t$ assuming the system starts with $m$ molecules. That is, $p_{n|m}$ is given by Eq. 22 with initial condition $p_{n}(0)=\delta_{nm}$. Eq. 23 becomes

[TABLE]

where the second step uses orthonormality, $\sum_{j}v_{nj}u_{jn^{\prime}}=\delta_{nn^{\prime}}$, and probability conservation, $u_{0n}=1$, to recognize that the $j=0$ term evaluates to $\bar{n}^{2}$. Inserting Eq. 25 into Eq. 9 and performing the integral (recalling that $\lambda_{j}<0$ for $j>0$), we obtain

[TABLE]

In matrix notation,

[TABLE]

where $\vec{n}$ is a row vector, $\vec{w}=mv_{m0}$ is a column vector, and neither the eigenvector matrices ${\bf V}$ and ${\bf U}$ nor the diagonal matrix $F_{jj^{\prime}}=-\delta_{jj^{\prime}}/\lambda_{j}$ contain the $j=0$ term. Numerically, we compute $\tau_{c}$ via Eq. 27 using a cutoff $N>n_{c}$ for the vectors and matrices.

The second method is to calculate $\tau_{c}$ from stochastic simulations Gillespie (1977) and the method of batch means Thompson (2010). The idea is to divide a simulation trajectory of length $T$ into batches of length $\tau_{b}$. In the limit $T\gg\tau_{b}\gg\tau_{c}$, the correlation time can be estimated by Thompson (2010)

[TABLE]

where $\sigma_{b}^{2}$ is the variance of the means of the batches.

In Fig. 5 we verify that the two methods converge to the same limit for sufficiently large $N$ or $T$, respectively. We find that the first method is more efficient until $n_{c}\sim 1000$, when numerically computing the eigenvectors for large $N>n_{c}$ becomes intractable.

Appendix B Experimental methods

The data in Fig. 6 [of which the smallest and largest doses are reproduced in Fig. 3(d)] were acquired at the same time and in a similar way as the data published in Vogel et al. (2016) and summarized in Erez et al. (arXiv:1703.04194). The difference is that, instead of only recording the data after steady state was reached, the time series was sampled by applying a chemical fixative to stop chemical reactions and preserve all biomolecular states. Specifically, we administered ice cold formaldehyde in PBS to each experimental well of a 96 well-v-bottom plate such that the final working dilution is 2%, and then transferred the cell-fixative solution to a new 96 well-v-bottom plate on ice. Cells were kept on ice for 10 minutes and then precipitated by centrifugation, resuspended in ice-cold 90% methanol, and placed in a $-20$ oC freezer until measurements were taken.