Cells in multicellular organisms have the ability to actively control their own death and engage in an organized self-destruction process known as programmed cell death, or apoptosis. This ability is necessary in many situations, ranging from embryonic development to the maintenance of mature organisms, and its deregulation is a hallmark of cancer [1] and many other diseases [2].

Apoptosis can be triggered by a variety of internal or external stimuli, such as specific death signals, lack of growth factors, DNA damage, and different types of stress [3]. These stimuli frequently converge with the metabolic state of the cell into the mitochondrial pathway, which promotes apoptosis through activation of the caspase cascade by releasing cytochrome c into the cytosol [3–5]. Low intracellular ATP levels favor the release of cytochrome c, resulting in high apoptotic rates [6, 7]. In this way, the absence of nutrients or signals that stimulate metabolism triggers apoptosis in a wide variety of cell types (Figure 1A).

The coupling of metabolism with apoptosis is remarkably important in cancer, especially because many anticancer drugs targeted against growth promoting pathways selectively induce apoptosis in cells with highly upregulated metabolism [8–10]. Extensive work has been done along these lines to decipher how different cellular components are wired to coordinate apoptotic responses [11, 12]. The connection of the single cell behavior with the collective dynamics of cell populations, however, has remained largely unexplored.

It is still not clear how the death-or-alive binary decision of single cells leads to the observed analog response of the cell population, where death increases continuously with decreasing ATP levels [13]. This connection is particularly challenging because classic population dynamics approaches [14] would require the death rate to peak at an intermediate ATP level to account for adaptation at low ATP levels, which paradoxically would imply that low ATP levels strongly favor survival. Here, the single-cell and population scales are connected through a quantitative mathematical model that identifies cellular variability, or noise, as the key element.

In general, cells have different levels of ATP, of proteins that promote or prevent cytochrome c release, and of the components of the caspase cascade. The most straightforward assumption is to consider that for each individual cell there is a specific value of the ATP level, an ATP level threshold, below which apoptosis is triggered and that intercellular variability in the cellular components makes this threshold noisy. Thus, a cell will die whenever the threshold is crossed, either because its ATP level decreases below its threshold or, vice versa, because its threshold increases over its ATP level.

In this paper, the focus is on death promoting stimuli that strongly affect metabolism, such as growth factor or nutrient withdrawal. Under these conditions, intracellular ATP levels change much faster than thresholds do and it is assumed that the noise in the threshold is quenched. Explicitly, the threshold is considered to be a time-independent random variable with precise statistical properties (Figure 1B).

The connection between the threshold distribution and cell death can be expressed in mathematical terms through the survival ATP-function, S(a), which is defined as the fraction of cells with a threshold below the average intracellular ATP level, denoted here by a. The initial threshold distribution is therefore given by $T(x)={\frac{dS(a)}{da}|}_{a=x}$, where x is the ATP threshold.

which uses the chain rule for derivatives: dS (a)/dt = (dS (a)/da) (da/dt). When the ATP level increases with time, within this model, the threshold would never be crossed and the apoptotic death rate would be zero.

A far-reaching implication of the threshold crossing assumption is that the death rate is not a function of just the ATP level but also of the ATP rate of change: da/dt. Thus, a precise testable prediction of the model is that changing the ATP level from decreasing (da/dt < 0) to increasing (da/dt > 0) will bring the apoptotic death rate to zero, irrespective of the value of the ATP level itself.

This model implements that the population dynamics response to energy depletion depends on both the energetic state of the cell and its dynamics. The complexity of cell death regulation [12] is collapsed into the form of the function r(a, da/dt), which can be inferred from experimental data, as shown below.

This approach can be applied straightforwardly to analyze experimental data for T-cell apoptosis upon growth factor withdrawal. Following loss of receptor engagement, T-cells rapidly downregulate the glucose transporter, glut1, which leads to lower intracellular ATP levels and results in apoptosis [15]. Experimental data is available for three levels of expression of the antiapoptotic protein Bcl-x_L[15]. For the three levels of expression of Bcl-x_L, the average ATP levels follow very similar temporal dynamics, thus indicating that Bcl-x_L does not affect cellular metabolism substantially.

The cellular metabolism experimental data (Figure 2A) can be incorporated into the model by considering that the ATP level, a(t), decays exponentially from its initial value, a(0), to a basal value, ATP*, with a characteristic time τ _ATP , which implies a net ATP consumption given by

In this case, autophagy provides the energy that allows cells to maintain a basal metabolic state in the absence of nutrient uptake [16].

where here x is the threshold, and μ and σ are two parameters related to the mean, $\overline{x}=e^{\mu +\sigma /2}$, and variance, $\overline{x^2}-{\overline{x}}^2=(e^{\sigma }-1)e^{2\mu +\sigma }$, of the distribution. The coefficient of variation is given by ${(\overline{x^2}-{\overline{x}}^2)}^{1/2}/\overline{x}={(e^{\sigma }-1)}^{1/2}$ and can be viewed as a measure of the cell-death noise. In general, μ and σ depend on the experimental conditions, cell type, and expression levels of Bcl-x_L.

The model defined by Equations (1)-(5) accurately accounts for the observed apoptotic behavior in response to growth factor withdrawal (Figure 2B) when the two parameters that characterize the threshold distribution are appropriately chosen for each cell type.

The response includes a lag phase, extensive death, and adaptation. All these features, in contrast to the intracellular ATP dynamics, strongly depend on Bcl-x_L expression. For the lowest expression, the cell type "Bcl-x_L" shows substantial death as soon as the ATP level starts to drop, with a very short lag phase (~ 5 hours). In this case, there is no adaptation and after 25 hours the whole cell population is dead. In the case of the cell type "Bcl-x_L 1E1", the lag phase lasts 20 hours approximately and is succeeded by a period of substantial death (until a time around 75 hours) before adaptation settles in with a small subpopulation of surviving cells. With the highest expression, the cell type "Bcl-x_L 4.1" is very resistant to energy depletion. It has virtually no cell death until 40 hours after growth factor withdrawal and a large fraction of cells are able to survive.

The initial distributions inferred from the experimental data (Figure 3), by choosing the parameters μ and σ to reproduce the observed behavior, indicate that, as Bcl-x_L expression increases, the thresholds shift on average toward lower values of the ATP level, consistent with the well-known antiapoptotic properties of Bcl-x_L. At the same time, the distributions increase their relative spread, resulting in higher cell-to-cell variability for higher Bcl-x_L expression (Figure 4).

It is important to emphasize that other types of threshold distributions, besides lognormal, could also be used in the model. However, if they are able to accurately reproduce the experimental data, they will look very similar to the lognormal distributions of Figure 3 in the range where intracellular ATP changes (from 3.54 to 0.32 fmol/cell). The reason is that in this system the shape of the threshold distribution is completely determined by the death rate. This result follows straightforwardly from integration of Eq. (1), which leads to the relationship ${\displaystyle {\int }_{x_1}^{x_2}(r/(da/dt))da}=\ln (S(x_1)/S(x_2))$. For instance, a Gaussian distribution would be able to provide relatively accurate results for the cell type "Bcl-x_L" but would fail to do so for the other two cell types. A key feature of the lognormal distribution is its ability to recapitulate the experimental observations for the three cell types.

The threshold distribution of live cells evolves in time because cells die when their threshold is higher than the ATP level (Figure 5 and Additional File 1). For monotonously decreasing ATP levels, this time-dependent behavior can be expressed explicitly in mathematical terms as T(x, t) = (dS (x)/dx) Θ (a(t)-x), where Θ(a(t) - x) is the Heaviside step function, which is 1 when a(t) - x > 0 and 0, otherwise. The resulting distribution is normalized to the fraction of surviving cells, ${\displaystyle {\int }_0^{\infty }T(x,t)}dx=S(a(t))$.

The results of the model indicate that the increased variability of cells that overexpress Bcl-x_L leads to highly heterogeneous responses, allowing a subpopulation of cells to escape death. Thus, within this model, adaptation observed in the late stages of the response to energy depletion is not the result of the apoptotic machinery adapting to low ATP levels but the result of the combined effects of high threshold variability with intracellular ATP reaching a steady state. If there were no noise, the whole population would either die or survive. When all the cells die, there is no adaptation. When all the cells survive, there is no response to the death stimulus and therefore there is no adaptation either. Thus, control of the population dynamics requires controlling not only the average values of the factors that control death but also their statistical properties.

The precise details of the transition from the single-cell to the cell-population dynamics can be elucidated by considering finite populations explicitly (Figure 6). Mathematically, the number of live cells is given by $N_{live}(t)={\displaystyle {\sum }_{i=1}^N\Theta (a(t)-x_i)}$, where x _i is the threshold of the cell i of a population of N initial (live) cells. The threshold x _i is a random variable with the statistical properties given by the distribution for each cell type. At the single cell level, for N = 1, there is no adaptation: either the cell dies at some random time or survives. For small populations, whether there is adaptation depends on the realization. There is no adaptation if all or none of the cells survive. In the limit of large N, the continuous population dynamics, N _live (t)/N ≈ S (a(t)), is recovered.

Recently, there has been an increasing interest on the role that molecular noise can play on the cellular behavior. Noise, in the form of random fluctuations in the number of molecules has been observed in many cellular processes, such as gene expression [19–22], protein abundance regulation [23–25], intracellular oscillations [26–28], and behavioral variability [29, 30]. The results presented here indicate that the oncogene Bcl-x_L exploits the inherent stochastic nature of molecular events to generate cell-death variability, which transforms a single-cell discrete event like death into a graded population response and leads to adaptation.

These two noise-induced features have far-reaching physiological consequences. On the one hand, a graded death response is essential for a smooth control of total cell numbers in multicellular organisms [13, 31, 32]. On the other hand, adaptation is a key mechanism that allows a subpopulation of cells to survive in the face of death-promoting stimuli, leading to the eventual failure of many anticancer therapies [33]. In the case of T-cell apoptosis, these two apparently disconnected features are in fact two sides of a single control mechanism based on noisy thresholds to trigger death.