Before detailing the theoretical and experimental results, we first present several biological questions for which cell lineage statistics could provide essential insights.

Growth of individual cells is heterogeneous in a cellular population even under constant environmental conditions (Stewart et al., 2005; Wakamoto et al., 2005; Wang et al., 2010; Hashimoto et al., 2016). Whether genetic or non-genetic, such growth heterogeneity inevitably enables selection within a cellular population. Growth heterogeneity can increase the rate of a population’s growth compared to the mean replication (division) rate of individual cells, known as ‘growth rate gain’ (Hashimoto et al., 2016). Since population growth rate is one of the critical quantities that determine long-term evolutionary success, it is interesting to ask to what extent growth heterogeneity contributes to population growth rate and how the contributions change depending on cellular phenotypes, genotypes (e.g. species), and environmental conditions. Answering this question may uncover strategies of each organism regarding how it exploits inherent stochasticity for population growth.

As we detail below, a measure of selection strength, $S_{\mathrm{KL}}^{(1)}[D]$, can quantify the growth rate gain from growth heterogeneity. Furthermore, we show that one can quantitatively decompose $S_{\mathrm{KL}}^{(1)}[D]$ into the contributions of distinct characteristics of growth heterogeneity, such as variance and skewness of fitness distributions. In this study, we apply the cell lineage statistics framework to single-cell lineage data and unravel how the growth rate gain changes across environments and organisms.

When a population of cells faces environmental changes, response of individual cells can be uniform and heterogeneous (Lambert et al., 2014; Julou et al., 2020). In one scenario, individual cells might respond to an environmental change uniformly and contribute to the future population nearly equally with respect to the number of descendants. In another scenario, only a tiny fraction of the cell population could respond to an environmental change, and the descendants of the responders might dominate the entire future population. In this case, the selection within a population is intense, and the nature of a population’s response exclusively depends on these rare cell lineages. Typically, the responses of real cell populations would fall between these two extremes; it is therefore critical to ask how strongly selection occurs within cellular populations in response to environmental changes to understand their response and adaptation strategies.

The framework enables such quantification by evaluating the selection strength $S_{\mathrm{KL}}^{(1)}[D]$ of responding cell populations. Importantly, quantifying the selection strength $S_{\mathrm{KL}}^{(1)}[D]$ requires only the information of division counts in cellular lineages. Hence, the selection strength is measurable even for complex processes where clarifying the transitions of environmental conditions around cells is technically challenging. We indeed analyze cellular populations of E. coli regrowing from an early or late stationary phase and characterize distinct levels of selection depending on the duration of stationary phase.

Since various traits of individual cells, such as expression levels of particular genes (Elowitz et al., 2002), are heterogeneous in cellular populations, it is natural to ask how strongly trait heterogeneity correlates with the fitness of individual cell lineages. Quantifying such correlations will allow us to understand which traits are under strong selection and potentially crucial for long-term evolution.

The cell lineage statistics framework quantifies relationships between traits and fitness using fitness landscapes $h(x)$. Additionally, the overall correlation between the heterogeneity of traits and that of fitness can be quantified by the relative selection strength $S_{\mathrm{rel}}[X]$. In this study, we measure $h(x)$ and $S_{\mathrm{rel}}[X]$ for a growth-regulation sigma factor in E. coli to unravel whether its continuum expression level heterogeneity is correlated with the fitness heterogeneity of single cell lineages.

Clarifying trait and fitness correlations based on individual-cell-based analyses is difficult when growth and traits fluctuate rapidly over time and when the traits affect growth with delays. In such circumstances, instantaneous correlations between traits and growth might not report their relations correctly. On the other hand, the cell-lineage-based analysis of this framework can take the whole dynamics of traits in cell lineages into account. For example, if we expect that absolute expression levels are important for fitness, the expression level averaged in each cell lineage can be employed as the lineage trait, and its fitness landscape and selection strength are measurable. If large fluctuations affect cell fates and contribute to diversification of cell lineage fitness within a cellular population (Purvis and Lahav, 2013), the variances of expression levels can be taken as lineage traits, and one can evaluate their fitness landscape $h(x)$ and relative selection strength $S_{\mathrm{rel}}[X]$. Therefore, the assumption of a cell lineage as a unit of selection can significantly extend the choice of traits, including time-dependent properties, and can provide insights into cellular dynamics that cannot be gained without the lineage-based formulation of fitness and selection.

First, we briefly review the analytical framework of cell lineage statistics introduced in Nozoe et al., 2017. This framework allows us to quantitatively infer fitness differences associated with distinct states of cellular lineage traits and selection within a growing cell population from empirical single-cell lineage tree data. Time-lapse single-cell measurements provide cellular growth and division information in the form of lineage trees (Figure 1, Stewart et al., 2005). We regard a lineage $\sigma$ as a cell history traceable back from a descendant cell at the final time point $t=\tau$ (Figure 1B). For the case of cellular growth shown in Figure 1A, 22 cell lineages exist in the trees.

Figure 1

Download asset Open asset

We assign two types of probability weight to cellular lineages. One is retrospective probability, in which we assign equal weight ${\displaystyle P_{\mathrm{r}\mathrm{s}}(\sigma ):=1/N_{\tau }}$ to all lineages, where $N_{\tau }$ is the number of cells at the final time point ${\displaystyle t=\tau .P_{\mathrm{r}\mathrm{s}}(\sigma )}$ represents the probability of selecting the history of a cell present at the endpoints of lineage trees. Another is chronological probability, in which we assign the weight ${\displaystyle P_{\mathrm{c}\mathrm{l}}(\sigma ):=2^{-D(\sigma )}/N_0}$ to the lineages, where $D(\sigma )$ is the number of cell divisions on lineage $\sigma$ and N₀ is the initial number of cells at ${\displaystyle t=\tau .P_{\mathrm{c}\mathrm{l}}(\sigma )}$ represents the probability of choosing lineage $\sigma$ descending the tree from one of the ancestor cells at $t=0$ and selecting one branch with the probability $1/2$ at every cell division. $P_{\mathrm{rs}}(\sigma )$ and $P_{\mathrm{cl}}(\sigma )$ can be different in general when the number of cell divisions are variable among the cell lineages, as shown in Figure 1B.

We define retrospective and chronological probabilities for a lineage trait $X$ as ${\displaystyle Q_{\mathrm{r}\mathrm{s}}(x):=\sum _{\sigma :X(\sigma )=x}{P}_{\mathrm{r}\mathrm{s}}(\sigma )}$ and ${\displaystyle Q_{\mathrm{c}\mathrm{l}}(x):=\sum _{\sigma :X(\sigma )=x}{P}_{\mathrm{c}\mathrm{l}}(\sigma )}$, where $X(\sigma )$ is the value of trait $X$ for lineage $\sigma$. Here, we regard any measurable quantity associated with cellular lineages as a lineage trait $X$. For example, time-averaged expression levels and production rates of a drug-resistance protein were analyzed as lineage traits in the experiments of Nozoe et al., 2017. Intuitively, $Q_{\mathrm{cl}}(x)$ and $Q_{\mathrm{rs}}(x)$ represent the probabilities of finding the lineage trait value $X=x$ before and after selection, respectively.

Using these retrospective and chronological distributions, we define the fitness landscape for lineage trait $X$ as

where ${\displaystyle \mathrm{\Lambda }:=\frac{1}{\tau }\ln \frac{N_{\tau }}{N_0}}$ is the population growth rate. This definition relates the relative difference of the retrospective probability from the chronological probability to fitness. $h(x)$ becomes greater than $\tau \mathrm{\Lambda }$ if the lineage trait state $X=x$ is overrepresented in the retrospective probability relative to chronological probability and vice versa. Furthermore, if none of the states of lineage trait $X$ are overrepresented nor underrepresented, $h(x)$ becomes constant across the states and equals $\tau \mathrm{\Lambda }$ for all $x$. The fitness landscape $h(x)$ thus represents fitness differences mapped on the lineage trait space of $X$ (see Figure 2 and Box 1).

Figure 2

Download asset Open asset

One can also define ‘selection strength’ using $Q_{\mathrm{rs}}(x)$ and $Q_{\mathrm{cl}}(x)$ as

where ${\displaystyle J[Q_{\mathrm{c}\mathrm{l}}(X),Q_{\mathrm{r}\mathrm{s}}(X)]:=\sum _x\left(Q_{\mathrm{c}\mathrm{l}}(x)-Q_{\mathrm{r}\mathrm{s}}(x)\right)\ln \frac{Q_{\mathrm{c}\mathrm{l}}(x)}{Q_{\mathrm{r}\mathrm{s}}(x)}}$ is Jeffreys divergence, and ${\displaystyle ⟨h(X){⟩}_{\mathrm{r}\mathrm{s}}:=\sum _{x}h(x)Q_{\mathrm{r}\mathrm{s}}(x)}$ and ${\displaystyle ⟨h(X){⟩}_{\mathrm{c}\mathrm{l}}:=\sum _{x}h(x)Q_{\mathrm{c}\mathrm{l}}(x)}$ are the retrospective and chronological mean fitness for lineage trait $X$. Jeffreys divergence measures dissimilarity between two probability distributions. Therefore, $S_{\mathrm{JF}}[X]$ measures dissimilarity between the chronological and retrospective distributions caused by selection. Notably, one can link this dissimilarity to the difference in the mean fitness, as shown in Equation 2. Since Jeffreys divergence is non-negative, the retrospective mean fitness (mean fitness after selection) is equal to or greater than the chronological mean fitness (mean fitness before selection).

This measure of selection strength quantifies how strongly differences in the states of lineage trait $X$ correlate with the differences in lineage fitness. Therefore, one can unravel which traits correlate with lineage fitness strongly by evaluating this for traits of interest.

Likewise, we can define two alternative selection strength measures:

where ${\displaystyle D_{\mathrm{K}\mathrm{L}}[Q_{\mathrm{c}\mathrm{l}}(X)||Q_{\mathrm{r}\mathrm{s}}(X)]:=\sum _x{Q}_{\mathrm{c}\mathrm{l}}(X)\ln \frac{Q_{\mathrm{c}\mathrm{l}}(X)}{Q_{\mathrm{r}\mathrm{s}}(X)}}$ and ${\displaystyle D_{\mathrm{K}\mathrm{L}}[Q_{\mathrm{r}\mathrm{s}}(X)||Q_{\mathrm{c}\mathrm{l}}(X)]:=\sum _x{Q}_{\mathrm{r}\mathrm{s}}(X)\ln \frac{Q_{\mathrm{r}\mathrm{s}}(X)}{Q_{\mathrm{c}\mathrm{l}}(X)}}$ are the Kullback-Leibler divergence of the two distributions. Note that $S_{\mathrm{KL}}^{(1)}[X]+S_{\mathrm{KL}}^{(2)}[X]=S_{\mathrm{JF}}[X]$.

These three types of selection strength measures share identical properties in common: they are always non-negative and report the overall correlations between trait states and fitness. We exclusively used $S_{\mathrm{JF}}[X]$ as the selection strength measure in our previous study (Nozoe et al., 2017). However, $S_{\mathrm{KL}}^{(1)}[X]$, $S_{\mathrm{KL}}^{(2)}[X]$, and their difference $S_{\mathrm{KL}}^{(2)}[X]-S_{\mathrm{KL}}^{(1)}[X]$ possess their own unique biological meanings, as we detail in Results. We indeed evaluate both $S_{\mathrm{KL}}^{(1)}[X]$ and $S_{\mathrm{KL}}^{(2)}[X]$ for the empirical lineage data of various organisms and use them to unravel distinct effects of selection on fitness variances. Such meanings and roles of the different selection strength measures are clarified in this study.

Importantly, division count $D$ is also a lineage trait, and its selection strength sets the maximum bound for the selection strength of any lineage trait irrespectively of a choice of the selection strength measures as discussed in Appendix 3. Therefore, the selection strength relative to that of $D$ is bounded between 0 and 1 and evaluates how strongly the heterogeneity of $X$ correlates with the division count heterogeneity in a given cellular population. This relative measure is useful when comparing relative strength of correlations between lineage traits and fitness across conditions. In this study, we define relative selection strength as

and use it in the analysis.

All of the quantities introduced above are measurable without relying on any growth models. Thus, this cell lineage statistics framework is applicable to a wide range of single-cell lineage data.

To quantify contributions of growth heterogeneity to population growth, we first rewrite the definition of the selection strength $S_{\mathrm{KL}}^{(1)}[X]$ (Equation 3) as follows:

This shows that population growth rates can be decomposed into chronological mean fitness and selection strength. In particular, when we take division count $D$ as a lineage trait, its fitness landscape is $\tilde{h}(d)=d\ln 2$ (Appendix 3), and ${⟨\tilde{h}(D)⟩}_{\mathrm{cl}}/\tau$ represents the mean division rate of cellular lineages without selection. $S_{\mathrm{KL}}^{(1)}[D]/\tau$, thus, represents growth rate gain caused by the growth heterogeneity among the cellular lineages in a cellular population. Therefore, evaluating $S_{\mathrm{KL}}^{(1)}[D]/\tau \mathrm{\Lambda }$ from single-cell lineage data provides information on the contribution of growth heterogeneity to population growth.

To further examine the connections between the disparate selection measures and elucidate their meaning, we define a function of a variable $\xi$ as

This is the cumulant generating function (cgf) of $h(x)$ with respect to the chronological distribution $Q_{\mathrm{cl}}$. We have $K_X(0)=0$, and from the definition of fitness landscape $h(x)$ (Equation 1), we find

When the radius of convergence of the Taylor expansion of $K_X(\xi )$ around $\xi =0$ is at least 1, $K_X(1)$ can be expressed as the series using the cumulants of a fitness landscape as

where ${\displaystyle {\kappa }_n^{(X)}:={\frac{d^n{K}_X(\xi )}{d{\xi }^n}|}_{\xi =0}}$ is the $n$-th order cumulant, satisfying ${\kappa }_1^{(X)}={⟨h(X)⟩}_{\mathrm{cl}}$, and ${\kappa }_2^{(X)}=\mathrm{Var}{[h(X)]}_{\mathrm{cl}}={⟨h{(X)}^2⟩}_{\mathrm{cl}}-{⟨h(X)⟩}_{\mathrm{cl}}^2$. Hence,

which shows that population growth rates can be expanded by the cumulants of a fitness landscape for any lineage trait $X$. Additionally, since ${\kappa }_1^{(X)}={⟨h(X)⟩}_{\mathrm{cl}}$, comparing (Equation 8) and (Equation 12) yields

Therefore, $S_{\mathrm{KL}}^{(1)}[X]$ represents the total contribution of second and higher order cumulants to population growth.

The cumulant expansion allows us to quantify the relative contributions of various statistical features of fitness distributions to population growth, such as mean, variance, and skewness. We define the cumulative contribution up to the $n$-th order cumulant as

and note that $W_n^{(X)}$ converges to 1 as $n\to \mathrm{\infty }$. In particular, $W_1^{(X)}=\frac{{⟨h(x)⟩}_{\mathrm{cl}}}{\tau \mathrm{\Lambda }}$ and $W_2^{(X)}=\frac{1}{\tau \mathrm{\Lambda }}\left({⟨h(X)⟩}_{\mathrm{cl}}+\frac{1}{2}\mathrm{Var}{[h(X)]}_{\mathrm{cl}}\right)$. We will indeed measure $W_n^{(D)}$ for various cellular species under steady and non-steady environments in the experimental sections below.

The function $K_X(\xi )$ defined in (Equation 9) is useful as it provides various forms of fitness and selection measures by simple algebraic calculation, as shown in Table 1. In general, evaluating $K_X(\xi )$ and its derivatives at $\xi =0$ and $\xi =1$ gives the information of chronological and retrospective statistics, respectively (Appendix 3). Therefore, $K_X(\xi )$ contains complete information on the fitness distributions in both chronological and retrospective statistics.

The difference between the two selection strength measures $S_{\mathrm{KL}}^{(1)}[X]$ and $S_{\mathrm{KL}}^{(2)}[X]$ is determined by the higher order cumulants by the relation

(Appendix 3). When fourth or higher order cumulants are negligible, the third-order fitness cumulant ${\kappa }_3^{(X)}$, that is the skewness of fitness distribution, determines which selection strength measure is greater.

The relations among the fitness and selection strength measures can be graphically depicted by plotting $K_X^{\prime }(\xi )$ in the interval $0\le \xi \le 1$ (Figure 3A). $S_{\mathrm{KL}}^{(1)}[X]$ corresponds to the area between $y={⟨h(X)⟩}_{\mathrm{cl}}$ and $y=K_X^{\prime }(\xi )$; and $S_{\mathrm{KL}}^{(2)}[X]$ corresponds to the area between $y=K_X^{\prime }(\xi )$ and $y={⟨h(X)⟩}_{\mathrm{rs}}$ (Figure 3A). Therefore, the skewness of fitness distribution primarily determines the convexity of $K_X^{\prime }(\xi )$ (Figure 3B–G).

Figure 3

Download asset Open asset

The difference between the two selection strength measures can reveal the effect of selection on fitness variances. The slope of the tangent lines to $K_X^{\prime }(\xi )$ at $\xi =0$ and 1 corresponds to the chronological and retrospective fitness variances, respectively (Table 1). Therefore, when $K_X^{\prime }(\xi )$ is convex upward in the interval $0\le \xi \le 1$ (${\displaystyle {\kappa }_3^{(X)}<0}$, i.e., ${\displaystyle S_{\mathrm{K}\mathrm{L}}^{(1)}[X]>S_{\mathrm{K}\mathrm{L}}^{(2)}[X]}$, as in Figure 3D), the effect of selection is to decrease the lineage fitness variance in the retrospective distribution relative to the chronological distribution, whereas if $K_X^{\prime }(\xi )$ is convex downward (${\displaystyle {\kappa }_3^{(X)}>0}$, i.e., ${\displaystyle S_{\mathrm{K}\mathrm{L}}^{(1)}[X]<S_{\mathrm{K}\mathrm{L}}^{(2)}[X]}$, as in Figure 3B), selection increases the fitness variance. We indeed find cases of both kinds of behavior in the experimental lineage data, as will be seen below. Therefore, one can probe the effect of selection on fitness variances by comparing the two selection strength measures $S_{\mathrm{KL}}^{(1)}[X]$ and $S_{\mathrm{KL}}^{(2)}[X]$.

Significant differences between $S_{\mathrm{KL}}^{(1)}[X]$ and $S_{\mathrm{KL}}^{(2)}[X]$ indicate non-negligible contributions of higher-order cumulants. In such circumstances, the fitness distributions are far from Gaussian with significant skews or multiple peaks. Therefore, higher-order cumulants can also be used to probe the existence of sub-populations in cellular populations.

We mentioned above that the selection strength measure $S_{\mathrm{KL}}^{(1)}[D]$ represents growth rate gain caused by fitness heterogeneity. Likewise, another selection strength measure $S_{\mathrm{KL}}^{(2)}[D]$ represents a different consequence of fitness heterogeneity, that is, additional loss of growth rate under fitness perturbations.

From (Equation 1), and taking division count as a lineage trait, one can express population growth rate as

We now consider the response of population growth rate to perturbations that cause lineage fitness to change from $D(\sigma )\ln 2$ to $(1-ϵ)D(\sigma )\ln 2$, and rewrite the population growth rate as

We have $\mathrm{\Lambda }(0)=\mathrm{\Lambda }$, and note that $\mathrm{\Lambda }(ϵ)=\frac{1}{\tau }{K}_D(1-ϵ)$ from (Equation 9). Differentiating $\mathrm{\Lambda }(ϵ)$ with respect to $ϵ$, and evaluating at $ϵ=0$, we find

(see Appendix 3). This relation shows that the change of population growth rate for small $ϵ$ is proportional to the retrospective mean fitness of the unperturbed population. Since ${⟨\tilde{h}(D)⟩}_{\mathrm{rs}}=\tau \mathrm{\Lambda }+S_{\mathrm{KL}}^{(2)}[D]$ (Equation 4), the relative change of population growth rate is

Therefore, a population with higher selection strength will exhibit a greater change in population growth rate upon perturbation. The selection strength measure $S_{\mathrm{KL}}^{(2)}[D]$ represents additional loss of population growth rate due to division count heterogeneity before perturbation.

As we see below, one manifestation of $ϵ$ occurs via a cell removal operation. Consider the removal of a branch in the genealogical tree just after each cell division with the probability of $1-2^{-ϵ}$ (${\displaystyle ϵ>0}$) (Figure 4A). In this case, the probability that a cell remains in the population after a cell division is $2^{-ϵ}$, and the growth of cell lineages that originally divided $d$ times will be effectively reduced by the factor ${(2^{-ϵ})}^d$. Consequently, the number of cell lineages that reach the end time point will also be effectively reduced from $N_0\left({\sum }_d{2}^d{Q}_{\mathrm{cl}}(d)\right)$ to $N_0\left({\sum }_d{2}^{(1-ϵ)d}{Q}_{\mathrm{cl}}(d)\right)$. Therefore, the population growth rate under this branch removal operation is given by (Equation 17), and the relative change of population growth rate is

Figure 4 with 1 supplement see all

Download asset Open asset

We validated this relation by simulating population growth with and without the cell removal operation (Figure 4B–E and Figure 4—figure supplement 1). The result confirmed that the relative changes of population growth rates by the probabilistic removal of cells followed $-\left(1+\frac{S_{\mathrm{KL}}^{(2)}[D]}{\tau \mathrm{\Lambda }}\right)ϵ$ in all the conditions (Figure 4C–E). We also tested this relation for cell populations with positive mother-daughter correlations of division intervals, which are often found for eukaryotic cells (Nozoe and Kussell, 2020; Seita et al., 2021; Mosheiff et al., 2018; Kuchen et al., 2020). We confirmed that the response relation was valid irrespectively of the strength of mother-daughter correlations (Figure 4—figure supplement 1), which shows that the relation is general and independent of the specific dynamics of the cell division process.

In Appendices 1 and 2, we calculate the exact form of $K_D(\xi )$ for analytically-tractable models. We derive chronological and retrospective mean fitness, selection strength, and the cumulants of fitness landscapes from $K_D(\xi )$ to observe how the framework works. In particular, we show the analytical calculation for a cellular population in which cells divide with gamma-distributed uncorrelated interdivision times in Appendix 2 to understand the effect of inherent stochasticity on population growth. This analysis yields two conclusions: (1) Unlike the central limit theorem, the contribution of higher-order cumulants to population growth remains even in the long-term limit, and (2) the shape of the generation time distribution influences the cell population’s long-term growth rate by constantly introducing selection within the population. Therefore, the details of inherent stochasticity of interdivision times are essential for the long-term population growth rate.

Next, we apply this framework of cell lineage statistics to experimental single-cell lineage data of various organisms. The list includes bacterial cells (Escherichia coli and Mycobacterium smegmatis), unicellular eukaryotic cells (Schizosaccharomyces pombe), and mammalian cancer cells (L1210 mouse leukemia cells). This analysis aims to unravel whether the extent of growth rate gain from growth heterogeneity depends on the organisms and environments under constant growth conditions. As summarized in Tables 2 and 3, we used cellular lineage data newly obtained in this study as well as other existing datasets (Nozoe et al., 2017; Wakamoto et al., 2013; Nakaoka and Wakamoto, 2017; Seita et al., 2021). The E. coli and S. pombe data include several culture conditions to compare cumulants’ contributions to population growth across environments. The E. coli data were obtained using either agarose pad or the microchamber array microfluidic device, yielding genealogical tree information such as the one shown in Figure 1. The S. pombe and L1210 cell data were obtained with mother machine microfluidic devices (Wang et al., 2010), which provide isolated cell lineage information but discard tree information due to its cell exclusion scheme. We assumed that these isolated cell lineages would follow chronological statistics and evaluated chronological distributions and selection strength according to the method described in Materials and methods. All the data analyzed in this section were taken from cell populations growing at approximately constant rates.

First, we evaluated the first-order cumulants’ contributions $W_1^{(D)}=\frac{{\kappa }_1^{(D)}}{\tau \mathrm{\Lambda }}=\frac{{⟨\tilde{h}(D)⟩}_{\mathrm{cl}}}{\tau \mathrm{\Lambda }}$ (Equation 14), finding that ${\displaystyle W_1^{(D)}<1}$ for all the samples and conditions (Figure 5A). This result confirms that the chronological mean fitness cannot fully account for the population growth rates. This means that the division count heterogeneity present even in constant environments contributes to increasing the population growth rate. However, the extent of the contributions was different: $W_1^{(D)}$ for S. pombe was consistently closer to 1 than those for the other cell types except one condition (EMM, 34 °C), suggesting that S. pombe’s growth is less heterogeneous under most culture conditions.

Figure 5 with 2 supplements see all

Download asset Open asset

We next evaluated $W_2^{(D)}=\frac{{\kappa }_1^{(D)}+{\kappa }_2^{(D)}/2}{\tau \mathrm{\Lambda }}$, finding that $W_2^{(D)}\approx 1$ for most of the conditions (Figure 5A). This result indicates small contributions of the third or higher-order cumulants to population growth. Consistent with this result, $S_{\mathrm{KL}}^{(1)}[D]$ and $S_{\mathrm{KL}}^{(2)}[D]$ were almost identical in most conditions (Figure 5B). Note that $S_{\mathrm{KL}}^{(2)}[D]-S_{\mathrm{KL}}^{(1)}[D]$ depends only on the third or higher order cumulants (Equation 15). The chronological distributions $Q_{\mathrm{cl}}(D)$ of these samples were nearly symmetric in most cases; however, under the conditions where the deviations of $W_2^{(D)}$ from 1 are larger, such as S. pombe in EMM medium and L1210, the distributions were skewed slightly (Figure 5C–E and Figure 5—figure supplement 1). Such distribution skew was reflected in the convexity directions of $K_D^{\prime }(\xi )$-plots (Figure 5F–H and Figure 5—figure supplement 2). These results imply that cellular populations of S. pombe in EMM medium and of L1210 contain small subpopulations that follow distinct division statistics. In fact, it was previously demonstrated that the L1210 cell populations contain slow-cycling cell lineages that can survive for longer durations under exposure to an anticancer drug (Seita et al., 2021). Therefore, this analysis confirms that the differences in the two strength measures can be used for detecting subpopulations in cellular populations.

In S. pombe EMM medium conditions, $K_D^{\prime }(\xi )$ was convex downward in the interval $0\le \xi \le 1$ except for EMM 34°C (Figure 5F and Figure 5—figure supplement 2). Therefore, under certain conditions selection can increase fitness variance in the retrospective distributions relative to chronological distributions among cellular lineages.

We further applied the framework to the cell lineage data of E. coli populations regrowing from an early or late stationary phase. This analysis aims to uncover how strongly selection occurs upon environmental changes and whether the selection strength can differ under identical conditions depending on the conditions before regrowth. To conduct time-lapse observations of regrowing cell populations, we used a microfluidic device equipped with microchambers etched on a glass coverslip. We sampled E. coli cells either from an early or late stationary phase batch culture and enclosed the cells into the microchambers by a semipermeable membrane (Inoue et al., 2001; Hashimoto et al., 2016). We switched flowing media from stationary-phase conditioned medium to fresh medium at the start of time-lapse measurements and recorded the growth and division of individual cells (Figure 6A, see Materials and Methods).

Figure 6

Download asset Open asset

The growth curves reconstructed by counting the number of cells at each time point showed lags in regrowth (Figure 6B). The lag time was shorter for the populations from the early stationary phase. The lineage tree structures in the cell populations were markedly different between the conditions (Figure 6C and D). The tree structures were more uniform for the early stationary phase sample with multiple divisions in most cell lineages (Figure 6C), whereas those for the late stationary phase sample were more heterogeneous, with 90% of cells showing no divisions within the observation time (Figure 6D).

We analyzed these data and found $W_1^{(D)}=0.95\pm 0.02$ for the population from the early stationary phase and $W_1^{(D)}=0.27\pm 0.04$ for the population from the late stationary phase (Figure 6E). Therefore, the chronological mean fitness, ${⟨\tilde{h}(D)⟩}_{\mathrm{cl}}$, explains only 27% of the growth rate of the population regrowing from the late stationary phase. In other words, significantly strong selection occurred in the regrowth from the late stationary phase. We also found that $W_2^{(D)}\approx 1$ for the population from the early stationary phase, as observed for the E. coli populations growing at constant rates. In constrast, $W_2^{(D)}$ for the population from the late stationary phase was $0.61\pm 0.04$, and $W_n^{(D)}$ converged to 1 only after taking the cumulants up to approximately 10th-order into account (Figure 6E). This indicates a skew of the fitness distribution and validates the existence of subpopulations following distinct division statistics in the population from the late stationary phase in this time scale of regrowth (Figure 6F). Reflecting the extreme skew to the right of the chronological distributions $Q_{\mathrm{cl}}(D)$ (Figure 6F), $S_{\mathrm{KL}}^{(2)}[D]$ was significantly greater than $S_{\mathrm{KL}}^{(1)}[D]$ for the late stationary phase sample (Figure 6G).

These results indicate that the levels of selection in the regrowing processes strongly depend on the durations under stationary phase conditions. Therefore, the ability to quickly resume growth under favorable conditions is gradually lost in most cells in the stationary phase; only a fraction of cells in the population can contribute to the future cell population. However, we also remark that preserving such non-growing cell lineages can be beneficial when cell populations are exposed to harsh environments in unpredictable manners (Kussell and Leibler, 2005).

RpoS is a sigma factor that controls the transcription of a large set of genes (10% of the genome) in E. coli (Battesti et al., 2011). High RpoS expression usually correlates with growth suppression; RpoS is induced when cells enter stationary phases or encounter stress conditions, such as starvation, low pH, oxidative stress, high temperature, or osmotic stress. Elevated RpoS expression provokes the intracellular programs to shut down growth and resist the stress (Battesti et al., 2011). However, it remains poorly understood how the continuum heterogeneity of RpoS expression levels is linked to the lineage fitness and selection in exponentially growing cellular populations. We therefore applied the lineage statistics framework to the single-cell time-lapse data of an E. coli strain expressing an RpoS-mCherry fusion protein from the native chromosomal locus and green fluorescent protein (GFP) from a low copy plasmid.

We quantified the time-scaled fitness landscapes $h(X)/\tau$ and relative selection strength $S_{\mathrm{rel}}[X]$ (Equation 5) under three growth conditions, taking the time-averaged mean fluorescent intensity of RpoS-mCherry or GFP along each cell lineage (proxies of time-averaged intracellular concentrations) as $X$ (Figure 7). Since fluorescent intensity is a trait that takes continuous values, we binned the intensity values with the bin sizes around which selection strength values are relatively stable (Materials and methods). Furthermore, since the calculation of relative selection strength from empirical data always gives positive values, we compared the relative selection strength values with those calculated from the data in which the correspondences between division counts and trait values were randomized to confirm the confidence levels (Figure 7—figure supplement 1).

Figure 7 with 1 supplement see all

Download asset Open asset

The result shows that the fitness landscapes and selection strength of RpoS expression level differ significantly among the growth conditions (Figure 7). Under the glucose-37°C condition, the fitness landscapes of RpoS-mCherry and GFP expression were both decreasing functions (Figure 7A and B). Thus, high expression of RpoS-expression and GFP in an exponentially growing population are both linked with lower lineage fitness. However, while the fitness landscape of GFP expression were nearly constant and showed significant decrease of fitness only at high expression levels, the fitness landscape of RpoS-mCherry decreased steadily in the observed expression range (Figure 7A and B). Consequently, the relative selection strength for RpoS-mCherry was 2.6-fold larger than that for GFP (Figure 7C).

Under the glucose-30°C and glycerol-37°C conditions, the fitness landscapes for RpoS-mCherry level were also decreasing functions and close to each other but significantly downshifted from that for the glucose-37°C condition (Figure 7A). This result reveals that cells could have different fitness for the same expression levels of RpoS, depending on the growth conditions. The selection strength for RpoS-mCherry was larger than that for GFP under the glucose-37°C and glucose-30°C conditions (Figure 7C), which proves that the heterogeneity of RpoS expression in a population correlates with the lineage fitness more strongly than that of GFP under those conditions. On the other hand, the relative selection strength of RpoS-mCherry under the glycerol-37°C condition was the smallest among the three conditions and not significantly different from that of GFP (Figure 7C). This is due to the relatively flat fitness landscapes in the central ranges of the distributions $Q_{\mathrm{cl}}(x)$ (Figure 7A and B) and the smaller variations of $x$ in the population (Figure 7D and E). These results reveal that the continuum heterogeneity of RpoS expression level in a population does correlate with the lineage fitness, but its contribution to selection depends on growth conditions. In other words, the heterogeneity in the RpoS-mCherry expression levels can barely correlate with fitness heterogeneity under some conditions.

We also evaluated the contributions of fitness cumulants for RpoS-mCherry expression to the population growth rate. Under all the conditions, $W_1^{(X)}$ was lower than 1 (Figure 7F–H). Therefore, the contributions of the higher-order fitness cumulants are non-negligible. However, the deviation of $W_1^{(X)}$ from 1 for RpoS-mCherry under the glycerol-37°C condition was small (Figure 7H). Hence, in this growth condition, RpoS-mCherry expression barely correlated with fitness heterogeneity in the population.

Importantly, this analysis can simultaneously reveal the changes in fitness landscapes (Figure 7A) and chronological distributions (Figure 7D). Interestingly, the distributions of the RpoS-mCherry expression levels are close between the Glucose-37°C and the Glycerol-37°C conditions, but the fitness landscapes are close between the Glucose-30°C and the Glycerol-37°C conditions. These results imply that the distributions and the fitness landscapes may vary independently in different conditions. Therefore, cells can potentially modulate the selection strength in each environment either by changing the fitness landscape or by changing the distribution of expression levels.

We constructed and used a microchamber array for conducting single-cell time-lapse observation under controlled environmental conditions. A microchamber is a well etched on a glass coverslip. We used two types of microchamber array. One is an array of microchamber, whose dimension is 70 μm (w) × 55 μm (h) × 1 μm (d). This microchamber has a 21-μm×7-μm pillar for supporting the membrane in the middle. We used this microchamber array for the exponential-phase experiment of E. coli. Another is an array of microchamber, whose dimension is 40 μm (w) × 30 μm (h) × 1 μm (d). We used this type of microchamber array for the stationary-phase-regrowth experiment in Figure 6. We fabricated these microchamber arrays following similar procedures described in Hashimoto et al., 2016; Inoue et al., 2001.

The photomasks for the microchamber array were created by laser drawing (DDB-201-TW, Neoark) on mask blanks (CBL4006Du-AZP, CLEAN SURFACE TECHNOLOGY). The photoresist on mask blanks was developed in NMD-3 (Tokyo Ohka Kogyo). The uncovered chromium (Cr)-layer was removed in MPM-E30 (DNP Fine Chemicals), and the remaining photoresist was removed by acetone. Lastly, the slide was rinsed in MilliQ water and air-dried.

The microchamber array was created in glass coverslips by chemical etching. First, we coated a 1,000-angstrom Cr-layer on a clean coverslip (NEO Micro glass, No. 1., 24 mm × 60 mm, Matsunami) by evaporative deposition and AZP1350 (AZ Electronic Materials) by spin-coating on the Cr-layer. We transferred the photomask patterns using a mask aligner (MA-20, Mikasa). After developing the photoresist in NMD-3 and the Cr-layer in MPM-E30, the coverslip was soaked in buffered hydrofluoric acid solution (110-BHF, Morita Kagaku Kogyo) for 14 minutes 20 seconds at 23°C for glass etching. The etching reaction was stopped by soaking the coverslip in milliQ water. The remaining photoresist and the Cr-layer were removed by acetone and MPM-E30, respectively.

We used a polydimethylsiloxane (PDMS) pad to flow culture medium and control the environmental conditions around the cells in the microchamber array. The PDMS pad was designed to have a square bubble-trap groove, which prevents interference with bright-field microscopic imaging by air bubbles in flowing media.

To create a mold for the bubble-trap groove, we spin-coated SU-8 3050 (Kayaku Advanced Materials) on a silicon wafer (ID 447, $\varphi$ = 76.2 mm, University Wafer) and baked it at 95°C for 2 hr on a hot plate. The SU-8 layer was exposed to UV light on a mask aligner using a photomask and postbaked at 95°C for 2 hr. After cooled down to room temperature, the SU-8 photoresist was developed in the SU-8 developer (Kayaku Advanced Materials) and rinsed with isopropanol (Wako).

Part A and Part B of PDMS resin (SYLGARD 184 Silicone Elastomer Kit, DOW SILICONES) were mixed at 10:1 and poured onto the SU-8 mold. The air bubbles were removed under a decreased pressure for 30 min. The PDMS was cured at 65°C for 1 hour, and 20 mm × 20 mm square PDMS pad was cut out using a blade. We punched out two holes ($\varphi$ = 2 mm) in the PDMS pad for the inlet and outlet, and 10-cm silicone tubes (SR-1554, Tigers Polymer Corp., outer $\varphi$ = 2 mm, inner $\varphi$ = 1 mm) were inserted into the holes. The tubes were fixed to the holes by gluing a small amount of PDMS around the tubes at the holes. This PDMS pad was washed in isopropanol by sonication and autoclaved for the single-cell measurements.

We washed the microfabricated coverslips by sonication in contaminon (Wako), ethanol (Wako), and 0.1 M NaOH solution (Wako). The washed coverslips were rinsed in milliQ water by sonication and dried at 140°C for 30 min. The washed coverslip was soaked in 1% (v/v) 3-(2-aminoethylaminopropyl)trimethoxysilane solution (Shinetsu Kagaku Kogyo) for 30 min and incubated at 140°C for 30 min to create an amino group on the glass surface. The treated coverslip was washed in milliQ water for 15 min and dried at 140°C for 30 min. 1 mg NHS-LC-LC-Biotin (Funakoshi) was dissolved in 25 μl dimethyl sulfoxide and dispersed in 1 ml phosphate buffer (0.1 mM, pH8.0). A total of 200 μl of this biotin solution was placed on the coverslip and incubated at room temperature for 4 hr. The biotin solution was removed by soaking the coverslip in milliQ water.

We prepared a streptavidin-decorated cellulose membrane to enclose cells in the microchamber array while retaining a flexible environmental control. First, a 3 cm × 3 cm square cellulose membrane (Spectra/Por7 Pre-treated RC Tubing MWCO:25kD) was cut out and washed in milliQ water for 10 min. The membrane was incubated in a 0.1 M NaIO₄ solution with gentle shaking for 4 hr at 25°C. After the wash in milliQ water, the treated membrane was incubated in a 500-μl solution of streptavidin hydrazide (Funakoshi) (10 μg/ml, dissolved in 0.1 mM phosphate buffer (pH7.0)) with gentle shaking for 14 hr at 25°C. The membrane was again washed in milliQ water and stored at 4°C.

We used two E. coli strains: MG1655 and MG1655 F3 rpoS-mcherry (MG1655 ΔfliCΔfimAΔflu rpoS-mcherry/pUA66-PrplS-gfp). MG1655 was used in the regrowth experiment from the stationary phases (Figure 6). MG1655 F3 rpoS-mcherry was used for analyzing the growth in constant environments (Figures 5 and 7). In MG1655 F3 rpoS-mcherry, the three genes, fliC, fimA, and flu, were deleted, and mcherry gene was inserted downstream of rpoS gene to express RpoS-mCherry translational fusion protein. This strain also expresses green fluorescent protein (GFP) from a low-copy plasmid, pUA66-PrplS-gfp, taken from a comprehensive library of fluorescent transcriptional reporters (Zaslaver et al., 2006).

We used MG1655 F3 rpoS-mcherry E. coli strain and cultured the cells in M9 minimal medium (Difco) supplemented with 1/2 MEM amino acids solution (SIGMA) and 0.2% (w/v) glucose or glycerol as a carbon source. We set the cultivation temperature either at 37°C or 30°C.

To prepare E. coli cells for single-cell observation, we first inoculated a glycerol stock into a 3-ml culture medium and incubated it with shaking overnight under the same conditions of culture medium and temperature as those used in the time-lapse measurement. 30 μl of the overnight culture was inoculated in a 3-ml fresh medium and incubated with shaking until the optical density at $\lambda$ = 600 nm reaches 0.1-0.3. This exponential-phase culture was diluted to OD₆₀₀ = 0.05, and 0.5 μl of the diluted cell suspension was spotted on the microchamber array on a biotin-decorated coverslip. A 5-mm × 5-mm streptavidin-decorated cellulose membrane was placed gently on the cell suspension on the coverslip, and an excess cell suspension was removed by a clean filter paper. A small piece of agar pad made with the culture medium and 1.5% (w/v) agar was placed on the cellulose membrane to maintain the culture conditions around the cells until tight streptavidin-biotin bonding was formed between the coverslip and the membrane. After 5-min incubation, the agar pad was removed, and the PDMS pad for medium perfusion was attached on the coverslip via a square-frame two-sided seal (Frame-Seal Incubation Chambers, Bio-rad). We immediately filled the device with the fresh medium and connected it to a syringe pump on the microscope stage.

We used E. coli MG1655 strain and cultured the cells in Luria-Bertani (LB) medium at 37°C. To prepare the cells for the time-lapse experiment, a glycerol stock of this strain was inoculated into a 2 ml LB medium and cultured with shaking for 15 hours. The cell culture was diluted in 50 ml fresh LB medium to OD₆₀₀ = 0.005 and again cultured with shaking as a pre-culture. For preparing the early-stationary-phase conditioned medium, 7 ml pre-culture cell suspension at 8 hr (OD₆₀₀ ≈ 4.3) was spun down at 2600 G for 12 min. The supernatant was filtered through a 0.22-μm filter. For preparing cells for time-lapse observation, a 10-μl pre-culture cell suspension at 8 hr was mixed with 240 μl early-stationary-phase conditioned medium. One μl of this diluted cell suspension was placed on the microchamber array on a biotin-decorated glass coverslip. A 5-mm × 5-mm streptavidin-decorated cellulose membrane was placed gently on the cell suspension on the coverslip, and an excess cell suspension was removed by a clean filter paper. A small piece of a conditioned medium agar pad made with 1.5% (w/v) agar was placed on a cellulose membrane to maintain the early stationary phase condition during the incubation. After 5-min incubation, the conditioned medium agar pad was removed, and the PDMS pad for medium perfusion was attached on the coverslip via a square-frame two-sided seal. We immediately filled the device with the conditioned medium and connected it to a syringe pump. We maintained the chamber filled with the conditioned medium until we started the time-lapse observation. The conditioned medium was flushed away immediately before starting the time-lapse measurement by flowing fresh LB medium. After flowing 2 ml fresh LB medium at 32 ml/hr, the flow rate was decreased and maintained at 2 ml/hr throughout the time-lapse measurement.

We followed the same procedures for the late stationary phase sample except that we sampled the cells and prepared the conditioned medium from a 24-hr pre-culture cell suspension (OD₆₀₀ ≈ 3.0).

We used Nikon Ti-E inverted microscope equipped with Plan Apo $\lambda$ 100× phase contrast objective (NA1.45), ORCA-R2 cooled CCD camera (Hamamatsu Photonics), Thermobox chamber (Tokai hit, TIZHB), and LED excitation light source (Thorlabs, DC2100). The microscope was controlled by Micromanager (Edelstein et al., 2014). In the exponential phase experiments, we monitored 25-30 microchambers in parallel in one measurement and acquired the phase-contrast, RpoS-mCherry fluorescence, and GFP fluorescence images from each position with a 3-min interval. We repeated the time-lapse measurement for each culture condition three times. In the regrowth experiment from the stationary phases, we monitored 150-250 microchambers in parallel with a 3-min interval and acquired only phase-contrast images.

We analyzed the time-lapse images by ImageJ (Schneider et al., 2012). We extracted the information of cell size (projected cell area), RpoS-mCherry fluorescence mean intensity, and GFP fluorescence mean intensity of individual cells along with division timings on each cell lineage for the exponential phase experiment. We extracted only division timings on each cellular lineage for the regrowth experiments from the stationary phases and used this information for further analysis.

To observe how the framework works, we show the exact form of $K_D(\xi )$ for a class of discrete probability distributions containing Poisson, binomial and negative binomial distributions. Let $\overline{D}$ and $\overline{D}\varphi$ denote the mean and the variance of $Q_{\mathrm{cl}}(D)$ respectively (i.e., $\varphi$ is the Fano factor of division counts). When $Q_{\mathrm{cl}}(D)$ is Poisson, binomial or negative binomial distributions, $\overline{D}$ and $\varphi$ uniquely determine the form of probability distribution: $\varphi =1$ for Poisson; ${\displaystyle \varphi <1}$ for binomial; and ${\displaystyle \varphi >1}$ for negative binomial (Appendix 1—figure 1A). Then, $K_D(\xi )$ for these distributions have a closed form

Appendix 1—figure 1

Download asset Open asset

(Appendix 3). We then immediately obtain

Since ${lim}_{\varphi \to 2}{K}_D\left(1\right)=\mathrm{\infty }$, ${\displaystyle 0<\varphi <2}$ is the range that the Fano factor of division counts can take within this scheme.

Using (Equation 35) allows us to calculate the cumulative contribution of cumulants of a fitness landscape ${\displaystyle W_n^{(D)}}$ (Equation 14). Plotting $W_n^{(D)}$ shows that the contribution of higher-order cumulants becomes significant when $\varphi$ is large (Appendix 1—figure 1B). Also, evaluating the values of the derivative of Equation 35 at $\xi =0$ and $\xi =1$, we have

Therefore, $S_{\mathrm{KL}}^{\left(1\right)}\left[D\right]/\tau \mathrm{\Lambda }$ and $S_{\mathrm{KL}}^{\left(2\right)}\left[D\right]/\tau \mathrm{\Lambda }$ depend only on the Fano factor $\varphi$. In particular, $S_{\mathrm{KL}}^{\left(1\right)}\left[D\right]=S_{\mathrm{KL}}^{\left(2\right)}\left[D\right]$ has 2 roots ${\displaystyle \varphi =0,{\varphi }_0(=0.5857)}$; ${\displaystyle S_{\mathrm{K}\mathrm{L}}^{\left(1\right)}\left[D\right]>S_{\mathrm{K}\mathrm{L}}^{\left(2\right)}\left[D\right]}$ if ${\displaystyle 0<\varphi <{\varphi }_0}$ and ${\displaystyle S_{\mathrm{K}\mathrm{L}}^{\left(1\right)}\left[D\right]<S_{\mathrm{K}\mathrm{L}}^{\left(2\right)}\left[D\right]}$ if ${\displaystyle {\varphi }_0<\varphi <2}$ (Appendix 1—figure 1C). Plotting $K_D^{\prime }(\xi )$ confirms that the covexity direction changes around ${\varphi }_0$ (Appendix 1—figure 1D). These analyses demonstrate how one can extract detailed information regarding selection in populations from $Q_{\mathrm{cl}}\left(D\right)$.

To understand how inherent stochasticity affect long-term population growth rate and selection, we consider a cellular population in which cells divide stochastically following a probability distribution of generation times (interdivision times).

Let $g\left(x\right)$ and $z$ denote the probability density function of generation time $x$ and the mean number of offsprings per generation, respectively. We assume that the generation time correlation between parent and offspring can be ignored; i.e., $g(x)$ gives the probability density that offspring’s generation time becomes $x$. The Malthusian parameter $\lambda$ is the real root of the so-called Euler-Lotka equation (Fisher, 1930):

We remark that (Equation 39) also holds for correlated generation times such as Markov models (Lebowitz and Rubinow, 1974) by reinterpreting $g(x)$ as the probability distribution of generation times of parent cells across a steadily growing population. In such cases, $g(x)$ depends on $z$, and we cannot treat $g(x)$ in (Equation 43) independent of $z=2^{\xi }$. Here, we ignore any transgenerational correlations in generation time to illustrate the effect of the variation in generation time on $K_D(\xi )$ and selection strength measures with simple calculations. For this purpose, we further choose gamma distributions as $g(x)$, i.e.,