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

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

Tommy A. Byrd Department of Physics and Astronomy, Purdue University, West Lafayette, Indiana 47907, USA    Amir Erez Department of Molecular Biology, Princeton University, Princeton, NJ 08544, USA    Robert M. Vogel IBM T. J. Watson Research Center, Yorktown Heights, New York 10598, USA    Curtis Peterson Department of Physics and Astronomy, Purdue University, West Lafayette, Indiana 47907, USA Department of Physics and School of Mathematical and Statistical Sciences, Arizona State University, Tempe, Arizona 85287    Michael Vennettilli Department of Physics and Astronomy, Purdue University, West Lafayette, Indiana 47907, USA    Grégoire Altan-Bonnet Immunodynamics Group, Cancer and Inflammation Program, National Cancer Institute, National Institutes of Health, Bethesda, Maryland 20814, USA    Andrew Mugler Department of Physics and Astronomy, Purdue University, West Lafayette, Indiana 47907, USA
###### 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.

thanks: These authors contributed equally.thanks: These authors contributed equally.

## 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 is the molecular species of interest, and the other species , , , etc. form a chemical bath for [Fig. 1(a)]. Whereas previously we considered only the steady state distribution of 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 of , (ii) the response time of , and (iii) the driving time 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 in which is produced from bath species and/or itself (feedback),

 jrX+Y+r⇌(jr+1)X+Y−r, (1)

where are stoichiometric integers. The probability of observing molecules of species in steady state according to Eq. 1 is

 pn=p0n!n∏j=1fj, (2)

where is set by normalization, and is a nonlinear feedback function governed by the reaction network. The inverse of Eq. 2,

 fn=npnpn−1, (3)

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

 m≡n∗−ncnc,h≡2(fnc−nc)−f′′′ncn3c,θ≡2(1−f′nc)−f′′′ncn2c, (4)

respectively, where is defined by , and are the maxima of . Qualitatively, sets the typical molecule number, drives the system to a unimodal () or bimodal () state, and biases the system to high () or low () molecule numbers. The critical point occurs at . The state variables , , and scale according to the exponents , , , and 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 is either produced spontaneously from bath species , or in a trimolecular reaction from two existing molecules and bath species ,

 A\raisebox0.0ptto0.0pt\raisebox2.0pt$\ext@arrow0359\arrowfill@\relbar\relbar⇀k−1k+1$\ext@arrow3095\arrowfill@↽\relbar\relbark−1k+1X,2X+B\raisebox0.0ptto0.0pt\raisebox2.0pt$\ext@arrow0359\arrowfill@\relbar\relbar⇀k−2k+2$\ext@arrow3095\arrowfill@↽\relbar\relbark−2k+23X. (5)

In this case the feedback function is

 fn=aK2+s(n−1)(n−2)(n−1)(n−2)+K2, (6)

where we have introduced the dimensionless quantities , , and in terms of the reaction rates and the numbers of and molecules. Given Eqs. 4 and 6, the effective thermodynamic variables , , and can be written in terms of , , and or vice versa Erez et al. (arXiv:1703.04194), with 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 diverges at the critical point according to Pathria and Beale (2011)

 τc|h=0 ∼|θ|−νz, (7) τc|θ=0 ∼|h|−νz/βδ, (8)

where we expect for mean-field dynamics Hohenberg and Halperin (1977); Kopietz et al. (2010). Here the autocorrelation time is defined as

 τc=1κ(0)∫∞0dt κ(t), (9)

where is the steady-state autocorrelation function, is the variance, and we have taken the start time to be without loss of generality because the system is in steady state.

To confirm the value of , we plot vs.  at (Eq. 8). We calculate 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 indeed diverges with , and that the location of the divergence approaches the expected value as the molecule number increases. We also see that the height of the peak increases with due to the rounding of the divergence Stephens et al. (2013). The inset of Fig. 2(b) plots this dependence: we see that at the critical point scales like for large (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 is sufficiently large, falls off with with the expected scaling exponent of . 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 () and the dynamic expectation for mean-field systems ().

### 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 in terms of the dynamics of the mean molecule number as

 τr=1Δ¯n(0)∫tmax0dt Δ¯n(t), (10)

where the quench occurs at , we define , and we ensure that . We compute from the time-dependent distribution using the stochastic simulations. Examples of 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 and . Specifically, scales identically with as it does with (Eqs. 7 and 8), which defines the Euclidean distance from the critical point as

 dc=[(θβδ)2+h2]1/2. (11)

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 ( or ) independently of the other.

To test whether the response time increases with proximity to the critical point, we must define initial values and for the environment before the quench, and a series of values and for the environment after the quench that are varying distances from the critical point . 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 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 on is shown in Fig. 3(c), and we see that indeed decreases as 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 and field can be calculated via Eq. 4 Erez et al. (arXiv:1703.04194). The values of and 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 . The general shape of the and 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,

 τr=1Δ¯I(0)∫tmax0dt Δ¯I(t), (12)

where and min. We calculate the distance from criticality using Eq. 11 as before, here using the experimental values of and . We see in Fig. 3(f) that the response time decreases with the distance from criticality , 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 comes directly from the experimental data, the distance from criticality 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 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 increases with the entropy , 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 with entropy . Calculating the entropy in Fig. 3(h) requires a conversion between intensity and molecule number , 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 , setting either and , or and , where and 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, , 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 , and the system begins to lag. Applying this condition to Eqs. 7 and 8, and using the above expressions for and , one obtains

 θ ∼τ−1/(νz+1)d, (13) h ∼τ−βδ/(νz+βδ)d, (14)

respectively. Because or near criticality in the mean-field Ising class, we have

 m ∼τ−β/(νz+1)d, (15) m ∼τ−β/(νz+βδ)d, (16)

respectively. Therefore, if the system is driven at different timescales , the Kibble-Zurek mechanism predicts that plots of the rescaled variables vs.  or vs.  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 for large , where sets the typical number of molecules in the system. Therefore, we interpret as a “system size,” and we correct for finite-size effects in the following way. Combining the relation with Eqs. 7 and 8, and Eqs. 13 and 14, we obtain

 nc ∼τ2νz/(νz+1)d, (17) nc ∼τ2νz/(νz+βδ)d, (18)

for the driving of or , respectively. We choose arbitrarily for a particular driving time , and when we choose a new , we scale 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 is driven from a positive to a negative value, the bifurcation response is lagging, occurring at a value less than the critical value (supercooling). Conversely, when is driven from a negative to a positive value, the convergence occurs at a value greater than (superheating). In both directions, the lag is larger when the driving is faster, corresponding to smaller values of (from yellow to dark brown). We see in Fig. 4(b) that similar effects occur for the driving of . 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 (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 numerically from the master equation for 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

 ˙→p=L→p. (19)

where is a tridiagonal matrix containing the birth and death propensities for . The eigenvectors of L satisfy

 L→vj =λj→vj, (20) →ujL =λj→uj, (21)

where the eigenvalues obey with only vanishing for the steady state, and because L is not Hermitian Walczak et al. (2009). Because Eq. 19 is linear in , the solution is

 pn(t)=∑jn′ujn′pn′(0)eλjtvnj (22)

for initial condition . Calling and , we write the autocorrelation function (see Eq. 9) as

 κ(t)=−¯n2+∑mnpmnmn=−¯n2+∑mnpn|mpmmn, (23)

where is the steady-state distribution, and is the dynamic solution at time assuming the system starts with molecules. That is, is given by Eq. 22 with initial condition . Eq. 23 becomes

 κ(t) =−¯n2+∑mnmnvm0∑jujmeλjtvnj (24) =∑mnmnvm0∞∑j=1ujmeλjtvnj, (25)

where the second step uses orthonormality, , and probability conservation, , to recognize that the term evaluates to . Inserting Eq. 25 into Eq. 9 and performing the integral (recalling that for ), we obtain

 τc=1σ2∑mnmnvm0∞∑j=1ujm(1−λj)vnj. (26)

In matrix notation,

 τc=σ−2→nVFU→w, (27)

where is a row vector, is a column vector, and neither the eigenvector matrices and nor the diagonal matrix contain the term. Numerically, we compute via Eq. 27 using a cutoff for the vectors and matrices.

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

 τc=τbσ2b2σ2, (28)

where 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 or , respectively. We find that the first method is more efficient until , when numerically computing the eigenvectors for large 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 C freezer until measurements were taken.

## References

• Strogatz (2018) S. H. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering (CRC Press, 2018).
• Scheffer et al. (2009) M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin, S. R. Carpenter, V. Dakos, H. Held, E. H. Van Nes, M. Rietkerk, and G. Sugihara, Nature 461, 53 (2009).
• Scheffer et al. (2012) M. Scheffer, S. R. Carpenter, T. M. Lenton, J. Bascompte, W. Brock, V. Dakos, J. Van de Koppel, I. A. Van de Leemput, S. A. Levin, E. H. Van Nes, et al., Science 338, 344 (2012).
• Veraart et al. (2012) A. J. Veraart, E. J. Faassen, V. Dakos, E. H. van Nes, M. Lürling, and M. Scheffer, Nature 481, 357 (2012).
• Dai et al. (2012) L. Dai, D. Vorselen, K. S. Korolev, and J. Gore, Science 336, 1175 (2012).
• Wang et al. (2012) R. Wang, J. A. Dearing, P. G. Langdon, E. Zhang, X. Yang, V. Dakos, and M. Scheffer, Nature 492, 419 (2012).
• Sha et al. (2003) W. Sha, J. Moore, K. Chen, A. D. Lassaletta, C.-S. Yi, J. J. Tyson, and J. C. Sible, Proceedings of the National Academy of Sciences 100, 975 (2003).
• Meisel et al. (2015) C. Meisel, A. Klaus, C. Kuehn, and D. Plenz, PLoS Computational Biology 11, e1004097 (2015).
• Hohenberg and Halperin (1977) P. C. Hohenberg and B. I. Halperin, Reviews of Modern Physics 49, 435 (1977).
• Erez et al. (arXiv:1703.04194) A. Erez, T. A. Byrd, R. M. Vogel, G. Altan-Bonnet, and A. Mugler (arXiv:1703.04194).
• Kibble (1976) T. W. Kibble, J. Phys. A 9, 1387 (1976).
• Zurek (1985) W. H. Zurek, Nature 317, 505 (1985).
• Schlögl (1972) F. Schlögl, Zeitschrift für Physik 253, 147 (1972).
• Dewel et al. (1977) G. Dewel, D. Walgraef, and P. Borckmans, Zeitschrift für Physik B Condensed Matter 28, 235 (1977).
• Nicolis and Malek-Mansour (1980) G. Nicolis and M. Malek-Mansour, Journal of Statistical Physics 22, 495 (1980).
• Brachet and Tirapegui (1981) M. Brachet and E. Tirapegui, Physics Letters A 81, 211 (1981).
• Grassberger (1982) P. Grassberger, Zeitschrift für Physik B Condensed Matter 47, 365 (1982).
• Prakash and Nicolis (1997) S. Prakash and G. Nicolis, Journal of Statistical Physics 86, 1289 (1997).
• Liu et al. (2007) D.-J. Liu, X. Guo, and J. W. Evans, Physical Review Letters 98, 050601 (2007).
• Vellela and Qian (2009) M. Vellela and H. Qian, Journal of the Royal Society Interface 6, 925 (2009).
• Pathria and Beale (2011) R. K. Pathria and P. D. Beale, Statistical mechanics (Academic Press, 2011).
• Kopietz et al. (2010) P. Kopietz, L. Bartosch, and F. Schütz, Introduction to the functional renormalization group, vol. 798 (Springer, 2010).
• Gillespie (1977) D. T. Gillespie, The journal of physical chemistry 81, 2340 (1977).
• Thompson (2010) M. B. Thompson, arXiv preprint arXiv:1011.0175 (2010).
• Stephens et al. (2013) G. J. Stephens, T. Mora, G. Tkačik, and W. Bialek, Physical review letters 110, 018701 (2013).
• Savitzky and Golay (1964) A. Savitzky and M. J. Golay, Analytical Chemistry 36, 1627 (1964).
• Vogel et al. (2016) R. M. Vogel, A. Erez, and G. Altan-Bonnet, Nature Communications 7, 12428 (2016).
• Altan-Bonnet and Germain (2005) G. Altan-Bonnet and R. N. Germain, PLoS Biology 3, e356 (2005).
• Chandran et al. (2012) A. Chandran, A. Erez, S. S. Gubser, and S. Sondhi, Physical Review B 86, 064304 (2012).
• Elowitz et al. (2002) M. B. Elowitz, A. J. Levine, E. D. Siggia, and P. S. Swain, Science 297, 1183 (2002).
• Skoge et al. (2011) M. Skoge, Y. Meir, and N. S. Wingreen, Physical review letters 107, 178101 (2011).
• Mora and Bialek (2011) T. Mora and W. Bialek, Journal of Statistical Physics 144, 268 (2011).
• Del Campo and Zurek (2014) A. Del Campo and W. H. Zurek, International Journal of Modern Physics A 29, 1430018 (2014).
• Deutschländer et al. (2015) S. Deutschländer, P. Dillmann, G. Maret, and P. Keim, Proceedings of the National Academy of Sciences p. 201500763 (2015).
• Deffner (2017) S. Deffner, Physical Review E 96, 052125 (2017).
• Friedman et al. (2006) N. Friedman, L. Cai, and X. S. Xie, Physical Review Letters 97, 168302 (2006).
• Mugler et al. (2009) A. Mugler, A. M. Walczak, and C. H. Wiggins, Physical Review E 80, 041921 (2009).
• Shahrezaei et al. (2008) V. Shahrezaei, J. F. Ollivier, and P. S. Swain, Molecular Systems Biology 4, 196 (2008).
• Horsthemke and Lefever (1984) W. Horsthemke and R. Lefever, Noise-induced transitions (Springer, 1984).
• Cotari et al. (2013) J. W. Cotari, G. Voisinne, O. E. Dar, V. Karabacak, and G. Altan-Bonnet, Science Signaling 6, ra17 (2013).
• Erez et al. (2018) A. Erez, R. Vogel, A. Mugler, A. Belmonte, and G. Altan-Bonnet, Cytometry Part A 93, 611 (2018).
• Das et al. (2009) J. Das, M. Ho, J. Zikherman, C. Govern, M. Yang, A. Weiss, A. K. Chakraborty, and J. P. Roose, Cell 136, 337 (2009).
• Prill et al. (2015) R. J. Prill, R. Vogel, G. A. Cecchi, G. Altan-Bonnet, and G. Stolovitzky, PLoS ONE 10, e0125777 (2015).
• Lee et al. (2008) J. A. Lee, J. Spidlen, K. Boyce, J. Cai, N. Crosbie, M. Dalphin, J. Furlong, M. Gasparetto, M. Goldberg, E. M. Goralczyk, et al., Cytometry Part A 73, 926 (2008).
• Walczak et al. (2009) A. M. Walczak, A. Mugler, and C. H. Wiggins, Proceedings of the National Academy of Sciences 106, 6529 (2009).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters