Universality of biochemical feedback and its application to immune cells

TL;DR

This paper establishes a universal framework linking biochemical feedback in immune cells to the Ising model, enabling data-driven insights into cellular responses and molecule numbers without detailed molecular knowledge.

Contribution

It introduces a mapping of biochemical feedback models to the Ising model, allowing extraction of thermodynamic parameters from biological data.

Findings

Successful application to mouse T cell fluorescence data

Revealed differences in drug response mechanisms

Enabled inference of molecule numbers from fluorescence intensity

Abstract

We map a class of well-mixed stochastic models of biochemical feedback in steady state to the mean-field Ising model near the critical point. The mapping provides an effective temperature, magnetic field, order parameter, and heat capacity that can be extracted from biological data without fitting or knowledge of the underlying molecular details. We demonstrate this procedure on fluorescence data from mouse T cells, which reveals distinctions between how the cells respond to different drugs. We also show that the heat capacity allows inference of absolute molecule number from fluorescence intensity. We explain this result in terms of the underlying fluctuations and demonstrate the generality of our work.

Figures

Tables

Table 1. Table 1: Drugs and dose ranges of experimental data Vogel et al. ( 2016 ) . Doses are spaced logarithmically. Fig.  3 uses PD325901 and Dasatinib. Fig.  4 uses all drugs.

Drug	Inhibits	Dose range (nM)	Shape in Fig. 4
PD325901	MEK	$0.09$$-$$1000$	Up triangle
AZD6244	MEK	$2.4$$-$$5000$	Down triangle
Trametinib	MEK	$0.5$$-$$1000$	Right triangle
Dasatinib	SRC	$0.09$$-$$1000$	Circle
Bosutinib	SRC	$0.5$$-$$1000$	Square
PP2	SRC	$24$$-$$50,$$000$	Diamond

Equations

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

$${\dot{p}}_n=b_{n-1}{p}_{n-1}+d_{n+1}{p}_{n+1}-(b_n+d_n)p_n,$$

$$p_n=p_0\prod _{j=1}^n\frac{b_{j-1}}{d_j}=\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},$$

$$\frac{C}{k_{\rm B}}=-(1+\theta)\sum_{n=0}^{\infty}p_{n}(1+\log p_{n})\Bigg{(}\psi_{n}-\sum_{j=0}^{\infty}p_{j}\psi_{j}\Bigg{)},$$

$$\varnothing \underset{k_2}{\overset{k_1}{\rightleftharpoons }}X,2X\underset{k_4}{\overset{k_3}{\rightleftharpoons }}D,$$

$$D\stackrel{k_5}{\to }D+A,A\stackrel{k_6}{\to }A+X,A\stackrel{k_7}{\to }\varnothing ,$$

$$D\stackrel{k_8}{\to }D+B,B+X\stackrel{k_9}{\to }B,B\stackrel{k_{10}}{\to }\varnothing .$$

$$0=\dot{n}/k_2=c_0-n_{\ast }+c_2{n}_{\ast }^2-c_3{n}_{\ast }^3,$$

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

$$n_c=\frac{3}{2}+\frac{1}{6}\sqrt{3(4K^2-1)}.$$

$$f_{n_c}$$

$$f_{n_c}^{\prime }$$

$$f_{n_c}^{\prime \prime \prime }$$

$$K^2$$

$$s$$

$$a$$

$$f_n$$

$$n_c$$

$$f_{n_c}$$

$$f_{n_c}^{\prime }$$

$$f_{n_c}^{\prime \prime \prime }$$

$$K$$

$$s$$

$$a$$

$$\ell =\log (I_{1}n),n=\frac{e^{\ell }}{I_1}.$$

$$q(\ell )=\frac{p_n}{d\ell /dn}=n{p}_n.$$

$$f_n=\frac{n{p}_n}{p_{n-1}}=(n-1)\frac{n{p}_n}{(n-1)p_{n-1}}=(n-1)\frac{q(\ell )}{q(\tilde{\ell })},$$

$$\tilde{\ell }=\log [I_1(n-1)]=\log (I_{1}n)+\log (1-ϵ)\approx \ell -ϵ.$$

$$f_n$$

$$\approx n(1-ϵ)\left(1+ϵ\frac{q^{\prime }}{q}\right)\approx n\left[1+ϵ\left(\frac{q^{\prime }}{q}-1\right)\right]$$

$$=n+\frac{q^{\prime }}{q}-1,$$

$$\varphi (\ell )=f_n-n=\frac{q^{\prime }}{q}-1=Q^{\prime }.$$

$${\partial }_n=I_1{e}^{-\ell }{\partial }_{\ell },$$

$$0$$

$$=I_1{e}^{-\ell }{\partial }_{\ell }(I_1{e}^{-\ell }{\partial }_{\ell }\varphi )=I_1^2{e}^{-\ell }\left(-e^{-\ell }{\varphi }^{\prime }+e^{-\ell }{\varphi }^{\prime \prime }\right)$$

$$=(I_1^2{e}^{-\ell }{)}^2\left({\varphi }^{\prime \prime }-{\varphi }^{\prime }\right).$$

$${\varphi }^{\prime \prime }({\ell }_c)={\varphi }^{\prime }({\ell }_c).$$

Full text

I Introduction

Positive feedback is ubiquitous in biochemical networks and can lead to a bifurcation from a monostable to a bistable cellular state Mitrophanov and Groisman (2008); Tkačik et al. (2012); Das et al. (2009); Vogel et al. (2016). Near the bifurcation point, the bistable state often reflects a choice between two accessible but opposing cell fates. For example, in T cells, the distribution of doubly phosphorylated ERK (ppERK) can be bimodal Vogel et al. (2016). ppERK is 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).

The bifurcation point is similar to an Ising-type critical point in physical systems such as fluids, magnets, and superconductors, where a disordered state transitions to one of two ordered states at a critical temperature Goldenfeld (1992). In fact, universality tells us that the two should not just be similar, they should be the same: because they are both bifurcating systems, both types of systems should exhibit the same critical scaling exponents and therefore belong to the same universality class Goldenfeld (1992). Although this powerful idea has allowed diverse physical phenomena to be united into specific behavioral classes, the application of universality to biological systems is still developing Mora and Bialek (2011); Munoz (2018); Salman et al. (2012); Brenner et al. (2015); Pal et al. (2014); Ridden et al. (2015); Qian et al. (2016); Hidalgo et al. (2014).

Biological tools such as flow cytometry, fluorescence microscopy, and RNA sequencing allow reliable experimental estimates of abundance distributions, inspiring researchers to seek to apply insights from statistical physics to biological data. In particular, recent studies have demonstrated that biological systems on many scales, from molecules Mora et al. (2010), to cells Kastner et al. (2015); Krotov et al. (2014); De Palo et al. (2017); Chen et al. (2012); Aguilar-Hidalgo et al. (2018); Wan and Goldstein (2018), to populations Bialek et al. (2014); Attanasi et al. (2014); Cavagna et al. (2017), exhibit signatures consistent with physical systems near a critical point. However, some of these studies have come under scrutiny because some of the signatures, particularly scaling laws, can arise far from or independent of a critical point Schwab et al. (2014); Touboul and Destexhe (2017); Newman (2005). Part of the problem is that the identification of appropriate scaling variables from data can be ambiguous, and one is often left looking for scaling relationships in an unguided way.

Typical approaches to the interpretation of abundance distributions include fitting to either detailed mechanistic models of the underlying reaction scheme, or to an effective description of the data such as a Gaussian or lognormal mixture model. The former approach is usually difficult to parameterize and difficult to generalize to other systems. The latter approach often suffers from numerical issues (the likelihood is unbounded and the expectation-maximization algorithm can lead to spurious solutions Biernacki et al. (2003)). Moreover, the vicinity of a bifurcation point is precisely where a mixture analysis is most likely to fail. In contrast, mapping to a statistical physics framework is expected to be universal, in the sense that the precise microscopic details of a broad range of biochemical models are unimportant near the bifurcation point, as they are coarse-grained rather than particular reaction parameters.

Here we provide a framework for mapping well-mixed stochastic models of biochemical feedback to the mean-field Ising model and apply it to published data on T cells. This allows us to extract effective thermodynamic quantities from experimental data without needing to fit to a parametric model of the system. This makes the theory applicable to a broad class of biological datasets without worrying about model selection or goodness-of-fit criteria. The theory provides insights on how T cells respond to drugs and reveals distinctions between one type of drug response and another. Furthermore, we find that one of the thermodynamic quantities (the heat capacity) provides a novel way to estimate absolute molecule number from fluorescence level in bifurcating systems. We demonstrate that our results can be extended to cases where feedback is indirect and discuss further extensions, including to spatiotemporal dynamics.

II Results

We consider a 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)]. The reactions of interest produce or degrade an $X$ molecule, can involve the bath species, and in principle are reversible. We allow for nonlinear feedback on $X$, meaning that the production of an $X$ molecule in a particular reaction might require a certain number of $X$ molecules as reactants. This leads to an arbitrary number of reactions of the form

[TABLE]

where in the $r$th reaction, $j_{r}$ are stoichiometric integers describing the nonlinearity, $k_{r}^{\pm}$ are the forward ($+$) and backward ($-$) reaction rates, and $Y_{r}^{\pm}$ represent bath species involved as reactants ($+$) or products ($-$). A simple and well-studied special case of Eq. 1 is Schlögl’s second model 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$ (i.e., $R=2$, $j_{1}=0$, $j_{2}=2$, $Y_{1}^{+}=A$, $Y_{2}^{+}=B$, and $Y_{1}^{-}=Y_{2}^{-}=\emptyset$).

We assume that molecules are well-mixed and that the numbers of bath molecules are constant. The latter assumption is equivalent to integrating out all species but $X$, such that the feedback on $X$ arises directly from $X$ itself (Eq. 1). However, in general the feedback will be indirect, with $X$ regulating dynamic species in the bath that in turn regulate $X$ (this is almost certainly the case in the T cells we study here). Therefore, we consider this more general case later in Section II.4 and show that the results discussed below remain unchanged.

The master equation for the probability of observing $n$ molecules of species $X$ according to Eq. 1 is

[TABLE]

where $b_{n}=\sum_{r=1}^{R}J_{rn}^{+}$ and $d_{n}=\sum_{r=1}^{R}J_{rn}^{-}$ are the total birth and death propensities, and $J_{rn}^{+}=k_{r}^{+}n_{r}^{+}n!/(n-j_{r})!$ and $J_{rn}^{-}=k_{r}^{-}n_{r}^{-}n!/(n-j_{r}-1)!$ are the forward and backward propensities of each reaction pair. Here $n_{r}^{\pm}$ are the numbers of molecules of the bath species involved in reaction $r$, and the factorials account for the number of ways that $X$ molecules can meet in a reaction. The steady state of Eq. 2 is Van Kampen (1992); Gardiner (1985)

[TABLE]

where $p_{0}^{-1}=\sum_{n=0}^{\infty}(1/n!)\prod_{j=1}^{n}f_{j}$ is set by normalization. In the second step of Eq. 3 we define an effective birth propensity $f_{n}\equiv nb_{n-1}/d_{n}$ corresponding to spontaneous death with propensity $n$ [Fig. 1(a)]. In general, $f_{n}$ is an arbitrary, nonlinear feedback function governed by the reaction network. For the Schlögl model, it is $f_{n}=[aK^{2}+s(n-1)(n-2)]/[(n-1)(n-2)+K^{2}]$, 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}^{-}$. As a ubiquitous example we also consider the Hill function $f_{n}=a+sn^{H}/(n^{H}+K^{H})$ with coefficient $H$. Importantly, the inverse of Eq. 3,

[TABLE]

allows calculation of the feedback function from the distribution Walczak et al. (2009), as utilized when analyzing the experimental data later in Section II.2.

The quantity $f_{n}-n$ determines the dynamic stability: there can be either one or two stable states $n_{*}$ [Fig. 1(b)], and the transition from a monostable to a bistable regime occurs at a bifurcation point [Fig. 1(c) inset]. These deterministic regimes correspond stochastically to unimodal and bimodal distributions $p_{n}$, respectively, with maxima at $n_{*}$, while the bifurcation point corresponds to a distribution that is flat on top [Fig. 1(c)].

II.1 Ising mapping and scaling exponents

To understand the scaling behavior near the bifurcation point, we expand the stability condition $f_{n_{*}}-n_{*}=0$ to third order around a point $n_{c}$ satisfying $f^{\prime\prime}_{n_{c}}=0$. This choice of $n_{c}$ eliminates the quadratic term in the dynamic forcing $f_{n}-n$, equivalent to eliminating the cubic term in an effective potential as in Ginzburg–Landau theory Kopietz et al. (2010). Defining the parameters

[TABLE]

the expansion $f_{n_{c}}+f^{\prime}_{n_{c}}(n_{*}-n_{c})+f^{\prime\prime\prime}_{n_{c}}(n_{*}-n_{c})^{3}/3!-n_{*}=0$ becomes $h-\theta m-m^{3}/3=0$. This expression is equivalent to the expansion of the Ising mean field equation $m=\tanh[(m+h)/(1+\theta)]$ for small magnetization $m$, where $\theta=(T-T_{c})/T_{c}$ is the reduced temperature, and $h$ is the dimensionless magnetic field Kopietz et al. (2010). Therefore, in our system we interpret $m$ as the order parameter, $\theta$ as an effective reduced temperature, and $h$ as an effective field. Explicit expressions for $n_{c}$, $\theta$, and $h$ in terms of the biochemical parameters and vice versa are given for the Schlögl and Hill models in Appendix A.

We see in Fig. 1(c) and (d) that $n_{c}$ determines where the distribution is centered, that $\theta$ drives the system to the unimodal ($\theta>0$) or bimodal ($\theta<0$) state, and that $h$ biases the system to high ($h>0$) or low ($h<0$) molecule numbers. Note that unlike in the Ising model, even when $h=0$ an asymmetry persists between the high and low states [see the purple distribution in Fig. 1(c)]. The reason is that in the master equation (Eq. 2), unlike in Ginzburg–Landau theory, fluctuations scale with molecule number, such that the high state is wider than the low state.

The equivalence between our system and the Ising mean-field equation near the critical point (Eq. 5) implies that our system has the same scaling exponents $\beta=1/2$, $\gamma=1$, and $\delta=3$ as the Ising universality class in its mean-field limit Kopietz et al. (2010). For completeness, we verify in Appendix B that these scalings are indeed obeyed by the Schlögl and Hill models.

However, Eq. 5 does not explicitly determine the value of the exponent $\alpha$. The reason is that, unlike $\beta$, $\gamma$, and $\delta$, the exponent $\alpha$ depends on the entire distribution $p_{n}$, not just the maxima. Specifically, $\alpha$ concerns the heat capacity, $C|_{h=0}\sim|\theta|^{-\alpha}$, which depends on the entropy $S$ and thus $p_{n}$. The equilibrium definition $C=T\partial_{T}S$ generalizes to a nonequilibrium system like ours when one uses the Shannon entropy $S=-k_{\rm B}\sum_{n}p_{n}\log p_{n}$ Mandal (2013). Since $T=(1+\theta)T_{c}$, we have $C=(1+\theta)\partial_{\theta}S$, or

[TABLE]

where $\psi_{n}\equiv(1/2)f^{\prime\prime\prime}_{n_{c}}n_{c}^{2}\sum_{j=1}^{n}(j-n_{c})/f_{j}$. Eq. 6 follows from performing the $\theta$ derivative using the expression in Eq. 3, the expansion below Eq. 5, and the definition of $\theta$ (Eq. 5). We see in Fig. 2(a) that when $h=0$, $C$ exhibits a minimum at $\theta^{*}$. We see in Fig. 2(b) that $\theta^{*}$ vanishes as the system size increases, $n_{c}\to\infty$. This implies that $C|_{h=0}\sim|\theta|^{0}$ to sub-quadratic order in $\theta$, or $\alpha=0$, again consistent with the Ising universality class in its mean-field limit. Interestingly, whereas $C$ is discontinuous in the mean-field Ising model Kopietz et al. (2010) and constant in the van der Waals model of a fluid Goldenfeld (1992), it is minimized here; nevertheless, in all cases $\alpha=0$. Note from Fig. 2(a) that $C$ is negative near $\theta=0$; negative heat capacity is a well-known feature of nonequilibrium steady states Zia et al. (2002); Boksenbojm et al. (2011); Bisquert (2005).

II.2 Application to immune cell data

To demonstrate the utility of our theory, we apply it to published data from T cells Vogel et al. (2016). In these experiments, chemotherapy drugs inhibit the enzymes MEK and SRC in the biochemical networks of the cells. The inhibition results in bimodal (low dose) or unimodal (high dose) distributions of ppERK abundance, which is measured as fluorescence intensity $I$ by flow cytometry. The distributions are shown for a range of drug doses in Fig. 3(a) and (b) (the insets show distributions of log intensity for clarity). Experimental details are given in the original publication Vogel et al. (2016) and are summarized in Appendix C, along with the drugs and dose amounts.

First, we compute the feedback function $f$ from each distribution using Eq. 4 (see Appendix D). Fig. 3(c) and (d) show the corresponding forcing functions [compare to Fig. 1(b)]. As expected, in each case we see that the forcing function transitions from two stable states to one stable state as the drug is applied.

Then, we compute $I_{c}$ (the analog of $n_{c}$ in units of fluorescence intensity), $\theta$, and $h$ from the feedback function using Eq. 5 (see Appendix D). These quantities are shown as a function of drug dose in Fig. 3(e)-(g). We see that the behavior is different depending on whether MEK inhibitor (MEKi) or SRC inhibitor (SRCi) is applied. Specifically, MEKi decreases $I_{c}$, increases $\theta$, and decreases $h$; whereas SRCi only decreases $h$, leaving the other quantities unchanged. Thus, the effective thermodynamic quantities can differentiate cellular responses to different perturbations, such as the application of different drugs.

Furthermore, the mapping provides an intuitive interpretation of the drug responses. MEKi causes a transition from a bimodal to a unimodal state in the expected way: by increasing the reduced temperature $\theta$ from a negative to a positive value [Fig. 3(f)]. In the process, $I_{c}$ decreases [Fig. 3(e)], meaning that the unimodal state is shifted to lower molecule number, near the lower mode of the bimodal state [Fig. 3(a) inset]. In contrast, SRCi causes a transition from a bimodal to a unimodal state in a different way: by decreasing the field while leaving $\theta$ and $I_{c}$ unchanged [Fig. 3(e)-(g)]. In essence, the distribution remains bimodal and unshifted, except that the field causes the high mode to diminish in weight [Fig. 3(b) inset]. Interestingly, the mean dose-response curves are similar for the two drugs Vogel et al. (2016), but our mapping elucidates precisely how the transitions are different at the distribution level. Related conclusions were drawn in Vogel et al. (2016), but those conclusions relied on fitting the distributions to a five-parameter Gaussian mixture model, which is expected to fail near the bifurcation point. Here we use only three parameters and no fitting, and we emerge with an intuitive interpretation in terms of thermodynamic quantities.

Finally, we note that for both drugs the effective field is negative at all doses [Fig. 3(g)]. The reason is that the fluorescence distributions have long tails (which is why they are often easier to visualize in log space); see Fig. 3(a) and (b). In the theory, a long tail is indistinguishable from a low-molecule-number bias in the peak, which corresponds to $h<0$. We address the possible origins and implications of the long tails in the Discussion (Section III).

II.3 Estimation of molecule number

We now apply the theory to compute the heat capacity from the T cell data. Specifically, we compute $C$ using Eq. 6 (see Appendix D) for all drugs and doses used in the experiments Vogel et al. (2016) (Appendix C). Unlike the other thermodynamic quantities, $C$ requires a conversion from fluorescence intensity to molecule number because it depends explicitly on the distribution $p_{n}$ (Eq. 6). Therefore we compute $C$ for various values of the conversion factor $I_{1}$, where $n=I/I_{1}$. The results are shown in Fig. 4. We see that irrespective of $I_{1}$ over four orders of magnitude, the data closest to $h=0$ (yellow) exhibit a global minimum in $C$ at $\theta=0$, as expected from Fig. 2(a). However, we also see that the depth of the minimum agrees with that of the theory only for the particular choice $I_{1}\approx 0.1$ [Fig. 4(c)].

To obtain a more precise estimate of $I_{1}$, we plot the sum of squared errors between the data and the theory as a function of $I_{1}$ in Fig. 4(e). We focus on the bifurcation region by considering only values of $\theta$ within $-\Delta\theta\leq\theta\leq\Delta\theta$, and we find that our results are not sensitive to the choice of $\Delta\theta$ [Fig. 4(f)]. This procedure (see the details in Appendix D) results in an estimate of $I_{1}=0.5$ $\pm$ $0.2$, as seen in Fig. 4(f). This value of $I_{1}$ corresponds to $\bar{n}_{*}=$ 170,000 $\pm$ 70,000 ppERK molecules in the high mode averaged across all cases with no inhibitor. It is possible to compare this value with previous measurements on these cells. In two separate experiments, it was estimated that there are approximately 100,000 Altan-Bonnet and Germain (2005) and 214,000 Hukelmann et al. (2016) ERK molecules per cell, and that only about 50% of these molecules are doubly phosphorylated during T cell receptor activation Altan-Bonnet and Germain (2005) (see Appendix D). These considerations give a range of roughly 50,000$-$107,000 ppERK molecules, which is consistent with our estimate of 170,000 $\pm$ 70,000. The agreement is especially notable given that T cell protein abundances generally span six orders of magnitude, from tens to tens of millions of molecules per cell Hukelmann et al. (2016).

Why does the heat capacity extract the conversion between fluorescence intensity and molecule number? As mentioned above, $\alpha$ is the only exponent that is a function of $p_{n}$ instead of just its maxima. This means that the plot of $C$ vs. $\theta$ contains information not only about means or modes, but also about fluctuations. The notion that fluctuation information is essential for converting from intensity to molecule number can be seen with a simpler example: a Poisson distribution. Here we would have $\sigma_{I}^{2}/\bar{I}^{2}=\sigma_{n}^{2}/\bar{n}^{2}=1/\bar{n}=I_{1}/\bar{I}$. From this relation it is clear that information about not only the mean ($\bar{I}$) but also the fluctuations ($\sigma_{I}^{2}$) in intensity is necessary and sufficient to infer the conversion factor $I_{1}$. In our case, the heat capacity is extracting similar information, but for a bifurcating system.

II.4 Generalization to indirect feedback

In the T cells, it is well known that ppERK does not apply feedback to its own activation directly, but rather indirectly via upstream components Vogel et al. (2016); Shin et al. (2009); Altan-Bonnet and Germain (2005). Therefore, we seek to determine the extent to which the above results are sensitive to our assumption in the theory that the feedback is direct. To this end, we construct a minimal extension of the model in Eq. 1 in which the feedback is indirect:

[TABLE]

Here $X$ is produced, is degraded, and reversibly dimerizes (first line); the dimer $D$ produces a species $A$ that produces $X$ and is degraded (second line); and the dimer also produces a species $B$ that degrades $X$ and is degraded (third line). Eq. 7 is an extension of Eq. 1 because there are multiple stochastic variables ($X$, $D$, $A$, and $B$), there are irreversible reactions, and $X$ feeds back on itself indirectly through $D$, $A$, and $B$ instead of directly.

The deterministic steady state of Eq. 7 is

[TABLE]

where $c_{0}\equiv k_{1}/k_{2}$, $c_{2}\equiv k_{3}k_{5}k_{6}/(k_{2}k_{4}k_{7})$, $c_{3}\equiv k_{3}k_{8}k_{9}/(k_{2}k_{4}k_{10})$, and the molecule numbers of $D$, $A$, and $B$ have been eliminated in favor of $n_{*}$ by setting their own time derivatives to zero. Because Eq. 8 is cubic in $n_{*}$, we see immediately that it has the same form as the expanded Ising mean field equation $h-\theta m-m^{3}/3=0$ (see Eq. 5). Specifically, defining $m=(n_{*}-n_{c})/n_{c}$ as in Eq. 5, the choice $n_{c}=c_{2}/(3c_{3})$ eliminates the term quadratic in $m$ and implies $\theta=3c_{3}/c_{2}^{2}-1$ and $h=9c_{0}c_{3}^{2}/c_{2}^{3}-3c_{3}/c_{2}^{2}+2/3$. It immediately follows that this model has the same exponents $\beta=1/2$, $\gamma=1$, and $\delta=3$ as the mean-field Ising universality class.

To test whether the heat capacity for this model exhibits the same features as that for the direct feedback model in Fig. 2(a), we compute the steady state marginal distribution $p_{n}$ using stochastic simulations Gillespie (1977) of Eq. 7. Specifically, we set $k_{3}/k_{4}=1/n_{c}$ and $k_{5}/k_{7}=k_{8}/k_{10}=1$ to ensure that the numbers of $D$, $A$, and $B$ molecules, respectively, are on the order of $n_{c}$. We then set $k_{4}/k_{2}=k_{7}/k_{2}=k_{10}/k_{2}=\rho$, where $\rho$ is a free parameter that determines whether the degradation timescales of $D$, $A$, and $B$, respectively, are faster ($\rho>1$) or slower ($\rho<1$) than that of $X$. These conditions, along with the definitions of $n_{c}$, $\theta$, and $h$ above, constitute nine equations for nine reaction rates, plus $k_{2}$ which sets the units of time. Solving these equations yields expressions for the rates in terms of $n_{c}$, $\theta$, $h$, and $\rho$ that we use in the simulations.

Fig. 5(a) shows the heat capacity $C$ as a function of $\theta$ for $h=0$, $n_{c}=100$, and $\rho=\{0.1,1,10\}$, where $C=(1+\theta)\partial_{\theta}S$ is computed from the entropy $S=-k_{\rm B}\sum_{n}p_{n}\log p_{n}$ by numerical derivative. We see that for all $\rho$ values, the curves exhibit a minimum at $\theta=0$, implying $\alpha=0$, and they rise more steeply for negative than for positive $\theta$ as in Fig. 2(a).

We then investigate whether Eq. 1 remains valid as a coarse-grained description of the extended model in Eq. 7. To answer this question, we infer values of $n_{c}$, $\theta$, $h$, and $C$ directly from the simulation data $p_{n}$ using the same protocol as for the experimental data. That is, we compute $f_{n}$ via Eq. 4, and then compute $\theta$, $h$, and $C$ from its derivatives at $n_{c}$ according to Eqs. 5 and 6, where $n_{c}$ satisfies $f^{\prime\prime}_{n_{c}}=0$. As with the experimental data (see Appendix D), derivatives are calculated using a Savitsky-Golay filter Savitzky and Golay (1964), although here we apply the filter directly to $f_{n}$ and perform the analysis directly in $n$ space, not log space.

Fig. 5(b) shows the result of this procedure for the inferred heat capacity $C$ as a function of the inferred $\theta$. We see that, as with the exact $C$ and $\theta$ [Fig. 5(a)], the data exhibit a minimum at $\theta=0$ and rise more steeply for negative than for positive $\theta$. Note that the values of $C$ and $\theta$ are different in (a) and (b), which is expected because the shape of $p_{n}$ is not quantitatively the same in the two models of Eqs. 1 and 7; nonetheless, the shape of the $C$ vs. $\theta$ curves remains the same. We have checked that the inferred values of $n_{c}$ and $h$ are distributed around their known values of $100$ and [math], respectively, and that the shape persists across a range of filter window sizes.

These results suggest that the main findings above are not sensitive to our assumption that feedback is direct, and therefore that we are justified in using Eq. 1 as a coarse-grained model to analyze the T cell data.

III Discussion

We have employed the fact that a feedback-induced bifurcation exhibits the scaling properties of the mean-field Ising universality class to provide a simple prescription for modeling and analyzing biological data. Contrary to existing mixture-model approaches, our method is most valuable near the bifurcation point, which is where biologically significant cell-fate decisions are expected to take place. Our approach provides the effective order parameter, reduced temperature, magnetic field, and heat capacity from experimental distributions without fitting or needing to know the molecular details. By applying the approach to T cell flow cytometry data, we discovered that these quantities discriminate between cellular responses in an intuitive, interpretable way, and that the heat capacity allows estimation of the molecule number from fluorescence intensity for a bifurcating system. By generalizing the theory to include indirect feedback, we demonstrated the capacity to model realistic signaling cascades where indirect feedback is common. Our approach should be applicable to other systems observed to undergo a pitchfork-like bifurcation and the associated unimodal-to-bimodal transition in abundance distributions, but not to systems which have an absorbing or extinction state, as they are expected to fall under a different universality class Ohtsuki and Keyes (1987); Grassberger and Sundermeyer (1978).

The theory assumes only birth-death reactions and neglects more complex mechanisms such as bursting Friedman et al. (2006); Mugler et al. (2009) or parameter fluctuations Shahrezaei et al. (2008); Horsthemke and Lefever (1984). These mechanisms are known to produce long tails and may be responsible for the long tails observed in the experimental data [Fig. 3(a) and (b)]. Cell-to-cell variability (CCV) may also contribute to the long tails, as it is known to be present in T cell populations Cotari et al. (2013). Our theory neglects CCV and instead assumes that the distribution of molecule numbers across the population is the same as that traced out by a single cell over time. Although CCV may play an important role, one generically expects the role of intrinsic fluctuations to be amplified near a critical point, and models that ignore CCV have been shown to be sufficient to explain both the bimodality Das et al. (2009) and variance properties Prill et al. (2015) of ppERK in T cells. Moreover, the fact that our theory provides an estimate of the molecule number that is consistent with other estimates suggests that intrinsic fluctuations play a large role. Distinguishing between intrinsic fluctuations and long-lived CCV is an important topic for future work.

Our work provides key tools that can be used for a broader exploration of biological systems. The approach is applicable to any experimental dataset that exhibits unimodal and bimodal abundance distributions, and could lead to a unified picture of diverse cell types and environmental perturbations in terms of effective thermodynamic quantities. At the same time, several extensions of our work are natural. For example, the dynamics of the theory could be probed to investigate the consequences of critical slowing down for driven or dynamically perturbed systems with feedback. Alternatively, the theory could be generalized to systems that are not well-mixed, such as intracellular compartments or communicating populations, to investigate space-dependent universal behavior and its biological implications.

IV Data availability

Data and code for all figures and the MIFlowCyt record are available at

https://github.com/AmirErez/UniversalImmune.

Acknowledgments

This work was supported by Human Frontier Science Program grant LT000123/2014 (Amir Erez), National Institutes of Health (NIH) grant R01 GM082938 (A.E.), Simons Foundation grant 376198 (T.A.B. and A.M.), and the Intramural Research Program of the NIH, Center for Cancer Research, National Cancer Institute.

Appendix A Mapping for Schlögl and Hill models

Here we provide the mapping from $n_{c}$, $\theta$, and $h$ to the biochemical parameters and vice versa for the Schlögl and Hill models. For the Schlögl model, the feedback function is

[TABLE]

The condition $f^{\prime\prime}_{n_{c}}=0$ is satisfied by

[TABLE]

The parameters $\theta$ and $h$ are given by Eq. 5, where

[TABLE]

These expressions are inverted to write the biochemical parameters $a$, $s$, and $K$ in terms of $n_{c}$, $\theta$, and $h$:

[TABLE]

where $x\equiv 2n_{c}-3$.

Similarly, for the Hill model we have

[TABLE]

In the Hill model, $H$ is an additional free parameter.

Appendix B Scaling exponents $\beta$, $\gamma$, and $\delta$

Here we verify that the stochastic Schlögl and Hill models have the scaling exponents $\beta$, $\gamma$, and $\delta$ of the mean-field Ising universality class. Specifically, we expect $m=\pm(-3\theta)^{\beta}$ for $h=0$ and $\theta<0$, with $\beta=1/2$; $\chi=\theta^{-\gamma}$ or $\chi=(-2\theta)^{-\gamma}$ for $\theta>0$ or $\theta<0$, respectively, with $\gamma=1$, where $\chi\equiv(\partial_{h}m)_{h=0}$ is the dimensionless susceptibility; and $m=(3h)^{1/\delta}$ for $\theta=0$, with $\delta=3$. Fig. 6 computes these quantities from the parameters and maxima of $p_{n}$ for the Schlögl and Hill models using the mapping in Eq. 5. We see that the scalings hold, as expected.

Appendix C Experimental methods

The experimental data analyzed in Fig. 3, along with a detailed description of the experimental methods, have been published previously Vogel et al. (2016). In this section we briefly summarize the experimental system and methods. The drugs and dose ranges used Figs. 3 and 4 are listed in Table 1.

The data investigate inhibition of the antigen-driven MAP kinase cascade in primary CD8+ mouse T cells. A natural way to stimulate T cells is to load a peptide (a fragment of an antigenic protein that the T cells are programmed to recognize) onto antigen-presenting cells. We achieve this by incubating RMA-S cells with antigen at $37$ oC. At the same time, we harvest the spleen and lymph nodes of a RAG2*-/-* OT1 mouse which has T cells specific only to the ovalbumin peptide with the amino acid sequence SIINFEKL. When we mix the OT1 T cells with the antigen-loaded RMA-S cells, we expose the OT1 T cells to their activating peptide. In response, the T cells activate their receptors through a SRC Family kinase (Lck). This triggers an enzymatic cascade, which in turn actives Ras-Raf-MEK-ERK leading to double phosphorylation of ERK, rendering it capable of communicating with the nucleus. By waiting for 10 minutes, the signaling reaches steady state and the distribution of the abundance of doubly phosphorylated ERK (ppERK) is the readout.

To measure the abundance of ppERK, we use fluorescence cytometry. Specifically, we introduce ppERK-targeted antibodies that are pre-conjugated with a fluorescent dye. Because antibodies selectively attach to their target molecule with negligible false-positives, the fluorescence intensity of the dye is proportional to the abundance of ppERK. To measure the intensity, approximately 30,000 cells per sample are passed one-by-one through a microfluidic device where they encounter a series of excitation lasers. Each cell yields one intensity value, and the histogram provides an estimate of the distribution of ppERK abundance across the population. We assume that the distribution across the population is a fair representation of the steady-state distribution of ppERK abundance of a single cell. This is reasonable (and is the accepted practice) since while the cells are alive and the experiment is taking place, they are in a dilute suspension (approximately 30,000 cells in 100 $\mu$L), not close enough together to influence each other.

Appendix D Experimental data analysis

We calculate the forcing functions and the effective thermodynamic quantities $I_{c}$, $\theta$, $h$, and $C$ from an experimental intensity distribution using the following procedure. First, we set $n=I/I_{1}$ to convert $p(I)$ to $p_{n}$, where the intensity of one molecule $I_{1}$ converts from intensity $I$ to molecule number $n$. We will see below that only $C$ will depend on the value of $I_{1}$.

Next, because the experimental distributions are long-tailed, we convert Eq. 5 to $\log I$ space for numerical stability. Here we provide the necessary conversions between functions of $n$ from the theory, and functions of $\ell\equiv\log I$ from the experiments, as the probability distributions over $n$ and $\ell$ do not have the same functional forms Erez et al. (2018). In what follows, prime denotes the derivative of a function with respect to its argument ($n$ for $f$; and $\ell$ for $q$, $Q$, and $\phi$). $\ell$ and $n$ are related as

[TABLE]

We denote the distribution of $\ell$ as $q(\ell)$. Approximating $n$ as continuous, probability conservation requires

[TABLE]

Using Eq. 26, the feedback function (Eq. 4) is

[TABLE]

where, using Eq. 25,

[TABLE]

The last steps define $\epsilon\equiv 1/n$ and assume that for most values of $n$ with appreciable probability we have $\epsilon\ll 1$. Therefore Eq. 27 becomes

[TABLE]

where we have kept to first order in $\epsilon$. Defining $\phi(\ell)\equiv f_{n}-n$, from Eq. 29 we have

[TABLE]

In the last step we define $Q(\ell)\equiv-\ell+\log q$ so that $\phi$ is computed as a total derivative, which we find more numerically stable. The $\phi(\ell)$ are the forcing functions plotted in Fig. 3(c) and (d).

The point $n_{c}$ is defined by $f^{\prime\prime}_{n_{c}}=0$. Eq. 25 implies

[TABLE]

such that the condition $f^{\prime\prime}_{n_{c}}=0$ becomes

[TABLE]

Therefore, we define a point $\ell_{c}$ by

[TABLE]

Numerically we enforce Eq. 33 by writing it as $0=\partial_{\ell}\left(\phi^{\prime}-\phi\right)$, and therefore

[TABLE]

Then

[TABLE]

from Eq. 25.

Derivatives of $f$ with respect to $n$ at $n_{c}$ are related in a straightforward way to derivatives of $\phi$ with respect to $\ell$ at $\ell_{c}$. First, the zeroth derivative is, by Eq. 30,

[TABLE]

where $n_{c}$ is defined in Eq. 35. Then, using Eq. 31, the first derivative is

[TABLE]

Finally, by a similar procedure, the third derivative is

[TABLE]

where the second step uses Eq. 33. Using Eqs. 35-D, $\theta$ and $h$ (Eq. 5) become

[TABLE]

Note that they do not depend on $I_{1}$.

To estimate the derivatives in Eqs. 39 and 40, we apply a Savitsky-Golay filter to the experimental $q(\ell)$ Savitzky and Golay (1964). Savitsky-Golay filtering replaces each data point with the value of a polynomial of order $J$ that is fit to the data within a window $W$ of the point. Since we require three derivatives of $\phi(\ell)$ (Eqs. 39 and 40), which depends on the first derivative of $q(\ell)$ (Eq. 30), we use the minimum value $J=4$. Thus, the procedure requires the adjustable parameter $W/L$, where $L$ is the number of $\log I$ bins. We find that $L=100$ and $W=25$ suffice [Figs. 3, 4, and 5(b)], and that results are robust to $W/L$.

The analysis is demonstrated for an example experimental distribution in Fig. 7. In summary, we:

plot $q(\ell)$ from the data using $L$ bins [Fig. 7(a), black]; 2. 2.

filter $q(\ell)$ using window $W$ [Fig. 7(a), red]; 3. 3.

compute $\phi(\ell)$ using Eq. 30 [Fig. 7(b)]; 4. 4.

compute $\ell_{c}$ using Eq. 34 [Fig. 7(c)]; 5. 5.

compute $I_{c}$, $\theta$, and $h$ from $\ell_{c}$, $\phi$, and its derivatives using Eqs. 35, 39, and 40; 6. 6.

compute $p_{n}$ from the data using $I_{1}$; and 7. 7.

compute $C/k_{\rm B}$ from $p_{n}$, $\theta$, $n_{c}$ (Eq. 35), $f^{\prime\prime\prime}_{n_{c}}$ (Eq. D), and $f_{n}$ (Eq. 29) using Eq. 6.

Fig. 7(d) shows that $I_{c}$ falls between the maxima as expected, and that $\theta$ and $h$ are negative corresponding to a distribution that is bimodal and skewed to the left, respectively.

To estimate the value of $I_{1}$, consider $\chi^{2}$, defined as

[TABLE]

where $N$ is the number of data points, $C_{i}/k_{\rm B}$ is the value of the heat capacity for each data point, $C(\theta_{i})/k_{\rm B}$ is the predicted value of the heat capacity at the location $\theta_{i}$ of that data point, and $\sigma_{i}^{2}$ is the variance for data point $i$. Under the simplifying assumption that $\sigma_{i}^{2}$ takes the same value $\sigma^{2}$ for all data points, we have $\chi^{2}=s/\sigma^{2}$, where $s$ is the sum of squared errors plotted in Fig. 4(e). As a function of $I_{1}$, $\chi^{2}$ should scale quadratically near its minimum,

[TABLE]

where the location and curvature of the minimum give the best estimate $\bar{I}_{1}$ and error in the estimate $\sigma_{I_{1}}$, respectively Bevington and Robinson (2002). In terms of $s$ we have

[TABLE]

where $s^{*}$ is the minimal value. The value $\sigma^{2}$ is, by definition, the average squared deviation of the data from the theory Bevington and Robinson (2002),

[TABLE]

here evaluated at the minimum $s^{*}$. Inserting this result into Eq. 43, we obtain

[TABLE]

We see that if $I_{1}$ deviates from $\bar{I}_{1}$ by $\sigma_{I_{1}}$, then $s$ is larger than its minimal value by a factor of $1+N^{-1}$. This criterion, illustrated by the black line in Fig. 4(e), is used to determine $\sigma_{I_{1}}$.

Fig. 4(e) is restricted to data whose $\theta$ values are less than or equal to $\Delta\theta=0.05$ in magnitude, of which there are $N=20$ points. As $\Delta\theta$ increases, $N$ increases, and the minimum of $s$ also becomes less sharp. These effects compensate, yielding an estimate of $I_{1}$ whose value and error are insensitive to $\Delta\theta$, as seen in Fig. 4(f). Averaged across $\Delta\theta$ values, we find $I_{1}=0.5$ and $\sigma_{I_{1}}=0.2$, as reported in the main text.

We compare our estimate of ppERK molecule number to two previous studies. In Altan-Bonnet and Germain (2005), it was estimated that there are 100,000 ERK molecules per cell (see Results in Altan-Bonnet and Germain (2005)). From Hukelmann et al. (2016), we estimate that there are 214,000 ERK molecules per cell. Specifically, from the Excel file associated with Fig. 1 in Hukelmann et al. (2016), we sum the mean number (column I) of ERK1 (also called MAPK3, row 2345) and ERK2 (also called MAPK1, row 874) to obtain 214,000 molecules to three significant digits. In Altan-Bonnet and Germain (2005), it was estimated that 50% of ERK molecules are doubly phosphorylated during T cell receptor activation (see caption of Fig. S2 in Altan-Bonnet and Germain (2005)).