A unified framework for measuring selection on cellular lineages and traits
Abstract
Intracellular states probed by gene expression profiles and metabolic activities are intrinsically noisy, causing phenotypic variations among cellular lineages. Understanding the adaptive and evolutionary roles of such variations requires clarifying their linkage to population growth rates. Extending a cell lineage statistics framework, here we show that a population’s growth rate can be expanded by the cumulants of a fitness landscape that characterize how fitness distributes in a population. The expansion enables quantifying the contribution of each cumulant, such as variance and skewness, to population growth. We introduce a function that contains all the essential information of cell lineage statistics, including mean lineage fitness and selection strength. We reveal a relation between fitness heterogeneity and population growth rate response to perturbation. We apply the framework to experimental cell lineage data from bacteria to mammalian cells, revealing distinct levels of growth rate gain from fitness heterogeneity across environments and organisms. Furthermore, third or higher order cumulants’ contributions are negligible under constant growth conditions but could be significant in regrowing processes from growtharrested conditions. We identify cellular populations in which selection leads to an increase of fitness variance among lineages in retrospective statistics compared to chronological statistics. The framework assumes no particular growth models or environmental conditions, and is thus applicable to various biological phenomena for which phenotypic heterogeneity and cellular proliferation are important.
Editor's evaluation
This manuscript presents a general statistical framework to infer selection on a quantitative trait, based on measurements of the values of this trait along related cell lineages. The manuscript provides both a detailed explanation of the mathematical underpinnings of the method and an illustration of its application to existing and new cell lineage datasets. This is a general framework and is not tailored to particular growth models or environmental conditions, making it applicable to broad examples of exponentially growing populations.
https://doi.org/10.7554/eLife.72299.sa0Introduction
Growth rates of cellular populations are physiological quantities directly linked to the fitness of cellular organisms. To understand the roles of biological processes and reactions within cells, including modulation of gene expression and metabolic states, one must characterize how they are eventually channeled into an increase or maintenance of population growth rates.
As documented by many singlecell studies, phenotypic states of individual cells in cellular populations are heterogeneous and often correlate with fitness variations among cellular lineages (Balázsi et al., 2011; Elowitz et al., 2002; Kelly and Rahn, 1932; Powell, 1956; Wakamoto et al., 2005; Wang et al., 2010; Cerulus et al., 2016; Susman et al., 2018). Fitness heterogeneity within a population causes a statistical bias on ancestral cells’ contributions to the number of descendants, which is broadly referred to as ‘selection’ (Leibler and Kussell, 2010). Such bias from growth heterogeneity makes the relations between cellular lineages and populations nontrivial. For example, an intriguing consequence of intrapopulation selection is a growth rate gain, a phenomenon that cell population’s growth rate becomes greater than the mean division rate of isolated singlecell lineages (Powell, 1956; Hashimoto et al., 2016; Rochman et al., 2018). Recent progress of singlecell measurements has enabled highthroughput acquisitions of cellular lineage trees and historical dynamics in each lineage (Stewart et al., 2005; Wang et al., 2010; Hashimoto et al., 2016). However, establishing the theory and method of cellular lineage statistics to quantify fitness differences among different phenotypic states and intrapopulation selection is still in progress (Nozoe et al., 2017; GarcíaGarcía et al., 2019; Levien et al., 2020; Genthon and Lacoste, 2020; Genthon and Lacoste, 2021).
Growth of cellular populations can be described using the ensemble of individual cells’ growth histories (Leibler and Kussell, 2010). A theoretical approach that regards a cell lineage (history) as a basic unit of analysis has offered illuminating insights into population dynamics. For example, it has provided the formula for untangling selection from responses (Leibler and Kussell, 2010), population response to agespecific changes in mortality and fecundity (Wakamoto et al., 2012), fluctuation relations of fitness (Kobayashi and Sughiyama, 2015; Genthon and Lacoste, 2020), and relations between cell size growth rate and population growth rate (Thomas, 2007; Lin and Amir, 2017).
Employing this cell historybased formulation of population dynamics, we have previously proposed a method of cellular lineage statistics that allows quantification of fitness landscapes and selection strength for any traits of cellular lineages (Nozoe et al., 2017). Here, we extend this statistical framework and show that population growth rates can be expanded by the cumulants that represent various properties of fitness distributions, such as variance and skewness, in a population. We apply the framework to experimental singlecell lineage data of bacteria, yeast, and mammalian cells to quantify their conditiondependent growth heterogeneity and its contribution to population growth rate. We also apply this framework to measuring the fitness landscapes for a growthregulating sigma factor in E. coli and identify the conditions where its continuum and nongenetic expression heterogeneity correlates with lineage fitness in cellular populations.
Examples of biological questions
Before detailing the theoretical and experimental results, we first present several biological questions for which cell lineage statistics could provide essential insights.
Growth rate gain
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 nongenetic, 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 longterm 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 singlecell lineage data and unravel how the growth rate gain changes across environments and organisms.
Selection in changing environments
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.
Correlations between cellular lineage traits and fitness
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 longterm 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 growthregulation 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 individualcellbased 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 celllineagebased 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 timedependent properties, and can provide insights into cellular dynamics that cannot be gained without the lineagebased formulation of fitness and selection.
Theoretical background
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 singlecell lineage tree data. Timelapse singlecell 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.
We assign two types of probability weight to cellular lineages. One is retrospective probability, in which we assign equal weight $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 $t=\tau .\phantom{\rule{thinmathspace}{0ex}}{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 $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_{0} is the initial number of cells at $t=\tau .\phantom{\rule{thinmathspace}{0ex}}{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 ${Q}_{\mathrm{r}\mathrm{s}}(x):=\sum _{\sigma :X(\sigma )=x}{P}_{\mathrm{r}\mathrm{s}}(\sigma )$ and ${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, timeaveraged expression levels and production rates of a drugresistance 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 $\mathrm{\Lambda}:=\frac{1}{\tau}\mathrm{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).
One can also define ‘selection strength’ using ${Q}_{\mathrm{rs}}(x)$ and ${Q}_{\mathrm{cl}}(x)$ as
where $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)\mathrm{ln}\frac{{Q}_{\mathrm{c}\mathrm{l}}(x)}{{Q}_{\mathrm{r}\mathrm{s}}(x)}$ is Jeffreys divergence, and $\u27e8h(X){\u27e9}_{\mathrm{r}\mathrm{s}}:=\sum _{x}h(x){Q}_{\mathrm{r}\mathrm{s}}(x)$ and $\u27e8h(X){\u27e9}_{\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 nonnegative, 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 $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)\mathrm{ln}\frac{{Q}_{\mathrm{c}\mathrm{l}}(X)}{{Q}_{\mathrm{r}\mathrm{s}}(X)$ and $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)\mathrm{ln}\frac{{Q}_{\mathrm{r}\mathrm{s}}(X)}{{Q}_{\mathrm{c}\mathrm{l}}(X)$ are the KullbackLeibler 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 nonnegative 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 singlecell lineage data.
A glossary of the terms
Here, we provide intuitive and illustrative explanations of the essential quantities in the cell lineage statistics and discuss their similarities and differences compared to the common usage in evolutionary biology.
Fitness: In evolutionary biology, fitness refers to the expected per capita contribution of individuals of a particular trait (usually a genotype) to the future population (Futuyma, 2010). For example, if a set of N_{0} individuals with trait $X$ produce N_{1} descendants on average in the future population, the fitness of this trait would be $N}_{1}/{N}_{0$. Since proliferation usually proceeds multiplicatively, the logarithm of fitness, In $({N}_{1}/{N}_{0})$, is also often referred to as ‘fitness’. Analogously, in our framework we define fitness for cell lineage traits as the expected contribution of lineages with a given trait value in the future population. For each cell lineage $\sigma$, the number of cell divisions occurring along the lineage, $D(\sigma )$, is used to estimate the expected contribution of each lineage to the future population.
Fitness landscape: In evolutionary biology, fitness landscapes are visual representations of relationships between reproductive abilities (fitness) and genotypes (Futuyma, 2010), where the height along the landscape corresponds to fitness. Since “genotype space” is vast and usually difficult to construct or visualize, fitness landscapes are often referred to as a metaphorical or conceptual tool for understanding complex evolutionary processes. For practical applications, however, fitness landscapes are often mapped on a low dimensional allele frequency space or a phenotypic space. Analogously, in our framework fitness landscapes are mapped on cell lineage trait spaces. However, they are different in that there is no assumption of genotypic differences underlying different trait states. Furthermore, the landscapes are directly measurable using cellular lineage trees and trait dynamics in each lineage.
For a cell lineage trait $X$, we define its fitness landscape to be a function $h(x)$ that reports the expected reproductive success of lineages having trait value $X=x$. Each lineage $\sigma$ having trait value $x$ contributes $2}^{D(\sigma )$ lineages to the future population, and by summing over lineages sharing the same trait value, we estimate the expected reproductive success of the trait and measure its fitness landscape. If differences in $X$ correlate with division count heterogeneity among cell lineages, $h(x)$ varies across the trait space of $X$; if differences in $X$ are uncorrelated with division count heterogeneity, $h(x)$ is constant over the entire space of $X$ (Figure 2).
Selection: The term selection refers to processes in which the frequencies of individuals with different traits change due to differences in their fitness (Futuyma, 2010). In evolutionary biology, selection is usually assessed based on changes in the distribution of traits between two points in times, which requires an accurate measure of fitness and a model to determine whether the observed changes were the result of trait fitness differences. In our cell lineage statistics framework, we measure selection by determining whether the observed distribution of lineage traits (i.e. the retrospective distribution) differs from the distribution expected in the absence of fitness differences (i.e. the chronological distribution). The key advantage that lineagebased analysis provides is the ability to construct explicitly the chronological distribution, which is the natural ‘null hypothesis’ against which selection can be tested in a modelindependent manner.
Selection strength: $S[X]$ (i.e., ${S}_{\mathrm{J}\mathrm{F}}[X]$, ${S}_{\mathrm{K}\mathrm{L}}^{(1)}[X]$, or ${S}_{\mathrm{K}\mathrm{L}}^{(2)}[X]$) is a quantitative measure that reports how strongly differences in the states of cell lineage trait $X$ are correlated with cell lineage fitness, taking the distributions of $X$ into account. The selection strength in our framework is measured by differences in the fitness measures or by differences between chronological and retrospective distributions (Equations 2–4). One can prove that these different definitions are mathematically equivalent.
The three situations depicted in Figure 2 would help us to gain an intuitive understanding of the properties and meanings of selection strength. When $X$ is correlated with fitness, a fitness landscape $h(x)$ becomes nonuniform, as mentioned above. When the states of lineage trait $X$ are heterogeneous and distributed widely within a population, the selection causes a significant difference between chronological and retrospective distributions due to the biased representation of trait states by selection. Therefore, the selection strength becomes large ($S[X]>0$, Figure 2A). In the second situation, $h(x)$ is again nonuniform, but the distribution of $x$ is narrow. In this case, there is almost no effective trait heterogeneity in the population on which selection can act. Consequently, the overall extent of selection becomes small, i.e., selection strength becomes small ($S[X]\approx 0$, Figure 2B). Finally, when $h(x)$ is uniform over the observed state of $x$, selection can neither overrepresent nor underrepresent any states, no matter how the trait $x$ distributes in a population. Therefore, the chronological and retrospective distributions become identical, and the selection strength becomes zero ($S[X]=0$, Figure 2C).
These examples show that $S[X]$ can gauge to what extent selection acts on a lineage trait $X$, considering both shapes of fitness landscapes and distributions of lineage traits in a population. Therefore, if $X$ is a trait of interest, quantifying $S[X]$ or the relative strength of selection $S[X]/S[D]$ determines how strongly the heterogeneity of $X$ is correlated with fitness differences of cell lineages.
In evolutionary biology, various measures are used to quantify how strongly selection acts in a population of interest. For example, the ‘coefficient of selection’ measures a relative difference in fitness of each genotype from that of the fittest genotype (Futuyma, 2010). This measure is useful when considering the selection against a particular reference genotype. The overall intensity of selection in a population can be quantified by changes in mean fitness before and after selection, variances of fitness before selection, changes in the mean of log fitness, and Jeffreys divergence between trait distributions before and after selection (Frank, 2012). Therefore, our definitions of selection strength follow the standard measures for the overall selection in evolutionary biology both conceptually and mathematically but are different in that the mean fitness and distributions of chronological and retrospective statistics are used.
Cumulants: In Results, we consider the contributions of the cumulants of a fitness landscape to population growth. The cumulants of a probability distribution are a set of quantities that characterize the distribution. For a discrete probability distribution $P(x)$, its cumulant generating function is defined as
and the $n$th order cumulant $\kappa}_{n$ is obtained by evaluating the $n$th order derivative of $K(\xi )$ at $\xi =0$, i.e.,
Notably, the first few cumulants correspond to important statistical quantities. The firstorder cumulant $\kappa}_{1$ corresponds to the mean $\u27e8X\u27e9:=E[X]=\sum _{x}xP(x)$, and the secondorder cumulant $\kappa}_{2$ corresponds to the variance $\mathrm{V}\mathrm{a}\mathrm{r}[X]:=E[{X}^{2}]E[X{]}^{2}=\sum _{x}{x}^{2}P(x){\left(\sum _{x}xP(x)\right)}^{2}$. The skewness of a distribution is usually defined as $E\left[{\left(\frac{XE[X]}{\sqrt{\mathrm{V}\mathrm{a}\mathrm{r}[X]}}\right)}^{3}\right]$, and this quantity can be expressed as $\kappa}_{3}/{\kappa}_{2}^{\frac{3}{2}$ using the thirdorder cumulant. Since $\kappa}_{2$ is positive, the sign of $\kappa}_{3$ determines the direction of the skewness: When ${\kappa}_{3}>0$, the distribution is skewed to the right with a long right tail; when ${\kappa}_{3}<0$, the distribution is skewed to the left with a long left tail.
Results
Growth rate gain and cumulant expansion of population growth rate
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 $\stackrel{~}{h}(d)=d\mathrm{ln}2$ (Appendix 3), and ${\u27e8\stackrel{~}{h}(D)\u27e9}_{\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 singlecell 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 $\kappa}_{n}^{(X)}:={\frac{{d}^{n}{K}_{X}(\xi )}{d{\xi}^{n}}}_{\xi =0$ is the $n$th order cumulant, satisfying ${\kappa}_{1}^{(X)}={\u27e8h(X)\u27e9}_{\mathrm{cl}}$, and ${\kappa}_{2}^{(X)}=\mathrm{Var}{[h(X)]}_{\mathrm{cl}}={\u27e8h{(X)}^{2}\u27e9}_{\mathrm{cl}}{\u27e8h(X)\u27e9}_{\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)}={\u27e8h(X)\u27e9}_{\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{{\u27e8h(x)\u27e9}_{\mathrm{cl}}}{\tau \mathrm{\Lambda}}$ and ${W}_{2}^{(X)}=\frac{1}{\tau \mathrm{\Lambda}}\left({\u27e8h(X)\u27e9}_{\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 nonsteady 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.
Difference in the selection strength measures reveals the effect of selection on fitness variance
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 thirdorder 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={\u27e8h(X)\u27e9}_{\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={\u27e8h(X)\u27e9}_{\mathrm{rs}}$ (Figure 3A). Therefore, the skewness of fitness distribution primarily determines the convexity of ${K}_{X}^{\prime}(\xi )$ (Figure 3B–G).
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$ (${\kappa}_{3}^{(X)}<0$, i.e., ${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 (${\kappa}_{3}^{(X)}>0$, i.e., ${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 nonnegligible contributions of higherorder cumulants. In such circumstances, the fitness distributions are far from Gaussian with significant skews or multiple peaks. Therefore, higherorder cumulants can also be used to probe the existence of subpopulations in cellular populations.
Population growth rate under fitness perturbations
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 )\mathrm{ln}2$ to $(1\u03f5)D(\sigma )\mathrm{ln}2$, and rewrite the population growth rate as
We have $\mathrm{\Lambda}(0)=\mathrm{\Lambda}$, and note that $\mathrm{\Lambda}(\u03f5)=\frac{1}{\tau}{K}_{D}(1\u03f5)$ from (Equation 9). Differentiating $\mathrm{\Lambda}(\u03f5)$ with respect to $\u03f5$, and evaluating at $\u03f5=0$, we find
(see Appendix 3). This relation shows that the change of population growth rate for small $\u03f5$ is proportional to the retrospective mean fitness of the unperturbed population. Since ${\u27e8\stackrel{~}{h}(D)\u27e9}_{\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 $\u03f5$ 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}^{\u03f5}$ ($\u03f5>0$) (Figure 4A). In this case, the probability that a cell remains in the population after a cell division is ${2}^{\u03f5}$, and the growth of cell lineages that originally divided $d$ times will be effectively reduced by the factor ${({2}^{\u03f5})}^{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\u03f5)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
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)\u03f5$ in all the conditions (Figure 4C–E). We also tested this relation for cell populations with positive motherdaughter 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 motherdaughter correlations (Figure 4—figure supplement 1), which shows that the relation is general and independent of the specific dynamics of the cell division process.
Applications to models
In Appendices 1 and 2, we calculate the exact form of ${K}_{D}(\xi )$ for analyticallytractable 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 gammadistributed 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 higherorder cumulants to population growth remains even in the longterm limit, and (2) the shape of the generation time distribution influences the cell population’s longterm growth rate by constantly introducing selection within the population. Therefore, the details of inherent stochasticity of interdivision times are essential for the longterm population growth rate.
Experimental evaluation of contributions of growth heterogeneity to population growth
Next, we apply this framework of cell lineage statistics to experimental singlecell 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 firstorder cumulants’ contributions ${W}_{1}^{(D)}=\frac{{\kappa}_{1}^{(D)}}{\tau \mathrm{\Lambda}}=\frac{{\u27e8\stackrel{~}{h}(D)\u27e9}_{\mathrm{cl}}}{\tau \mathrm{\Lambda}}$ (Equation 14), finding that ${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.
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 higherorder 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 slowcycling 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.
The contributions of higher order cumulants become significant in the regrowth from a late stationary phase
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 timelapse 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 stationaryphase conditioned medium to fresh medium at the start of timelapse measurements and recorded the growth and division of individual cells (Figure 6A, see Materials and Methods).
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, ${\u27e8\stackrel{~}{h}(D)\u27e9}_{\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 10thorder 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 nongrowing cell lineages can be beneficial when cell populations are exposed to harsh environments in unpredictable manners (Kussell and Leibler, 2005).
Lineage statistics reveal conditiondependent fitness landscapes and selection strength for a growthregulating sigma factor
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 singlecell timelapse data of an E. coli strain expressing an RpoSmCherry fusion protein from the native chromosomal locus and green fluorescent protein (GFP) from a low copy plasmid.
We quantified the timescaled fitness landscapes $h(X)/\tau $ and relative selection strength ${S}_{\mathrm{rel}}[X]$ (Equation 5) under three growth conditions, taking the timeaveraged mean fluorescent intensity of RpoSmCherry or GFP along each cell lineage (proxies of timeaveraged 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).
The result shows that the fitness landscapes and selection strength of RpoS expression level differ significantly among the growth conditions (Figure 7). Under the glucose37°C condition, the fitness landscapes of RpoSmCherry and GFP expression were both decreasing functions (Figure 7A and B). Thus, high expression of RpoSexpression 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 RpoSmCherry decreased steadily in the observed expression range (Figure 7A and B). Consequently, the relative selection strength for RpoSmCherry was 2.6fold larger than that for GFP (Figure 7C).
Under the glucose30°C and glycerol37°C conditions, the fitness landscapes for RpoSmCherry level were also decreasing functions and close to each other but significantly downshifted from that for the glucose37°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 RpoSmCherry was larger than that for GFP under the glucose37°C and glucose30°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 RpoSmCherry under the glycerol37°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 RpoSmCherry expression levels can barely correlate with fitness heterogeneity under some conditions.
We also evaluated the contributions of fitness cumulants for RpoSmCherry 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 higherorder fitness cumulants are nonnegligible. However, the deviation of ${W}_{1}^{(X)}$ from 1 for RpoSmCherry under the glycerol37°C condition was small (Figure 7H). Hence, in this growth condition, RpoSmCherry 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 RpoSmCherry expression levels are close between the Glucose37°C and the Glycerol37°C conditions, but the fitness landscapes are close between the Glucose30°C and the Glycerol37°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.
Discussion
Growth and division of individual cells are intrinsically variable, which causes division count heterogeneity among cellular lineages in a population. Such heterogeneity is ubiquitous across prokaryotic and eukaryotic cells, and its statistical properties could depend on the mechanisms and regulations determining cell division timings. Notably, division count heterogeneity influences population growth rate and, consequently, a population’s survival and evolutionary success. Therefore, understanding what statistical features are produced among cellular lineages and how these features contribute to population growth is essential for unraveling each organism’s survival and evolutionary strategy.
This report presents a cell lineage statistics framework to uncover the linkage between fitness distributions and population growth rate. We reveal that a population’s growth rate can be expanded by the cumulants of a fitness landscape for any lineage trait. The cumulant expansion allows us to quantify the contribution of each fitness cumulant, such as variance and skewness, to population growth rate. Applying this framework to the experimental cell lineage data revealed the cumulants’ contributions to population growth for various organisms and environmental conditions. In particular, higherorder cumulants became significant in the regrowth of E. coli from a late stationary phase. We remark that the cumulant expansion of population growth rate is valid only when all the cumulants are finite and when the Taylor expansion of ${K}_{X}(\xi )$ around $\xi =0$ also converges at $\xi =1$. However, all the experimental data examined in this study exhibited stable convergence, including in the regrowth condition from the late stationary phase.
An advantage of this framework is its independence from any growth and division models. The mechanisms driving the growth and division of individual cells are diverse among organisms. For example, the properties of cellular growth and division, such as whether a cell’s size increases exponentially or linearly and whether cell size regulation follows sizer or adder models, could depend on cell types, organisms, and environmental conditions (Jun et al., 2018; Kohram et al., 2021). Therefore, any model assumptions restrict applicability and necessitate model validation before application. The model independence of the framework presented here comes from the definitions of two essential quantities: the chronological and retrospective probabilities. Quantifying these probabilities requires only the information of the numbers of cells at initial and end time points and of division counts on each cellular lineage. Consequently, this formalism can be applied even to nonstationary conditions without modifications. However, we also remark that this independence from the details other than cell lineage structures imposes a limitation on the framework because it cannot report any potential influences from factors such as heterogeneous environments around cells and nonquantified traits. Furthermore, the fitness landscape $h(x)$ and the relative selection strength ${S}_{\mathrm{rel}}[X]$ evaluate only the correlations between the trait and fitness, not causal relationships. However, causal traits should have large selection strength values, and this framework helps narrow down the candidates for essential traits. Most importantly, division statistics is the focal information that connects molecular details underlying cellular growth and division to population growth. Regulatory mechanisms can influence population growth only by modulating the division statistics in a cellular population.
Growth heterogeneity in a cellular population plays a critical role in its adaptation and survival against stressful conditions. In antibiotic persistence, bacterial cell populations often harbor small populations of nongrowing or slowgrowing cells which can survive under antibiotic exposures (Balaban et al., 2004). Such structures of growth heterogeneity can be investigated in a unified manner by the selection strength measures introduced here. For example, the differences in ${S}_{\mathrm{KL}}^{(1)}[D]/\tau \mathrm{\Lambda}$ among organisms can reveal the distinct levels of the overall growth heterogeneity of these organisms. Furthermore, the differences between ${S}_{\mathrm{KL}}^{(1)}[D]$ and ${S}_{\mathrm{KL}}^{(2)}[D]$ characterize the structure of growth heterogeneity: If ${S}_{\mathrm{K}\mathrm{L}}^{(1)}[D]>{S}_{\mathrm{K}\mathrm{L}}^{(2)}[D]$, the distribution of lineage fitness is skewed negatively, and the cell population harbors small subpopulations of slowgrowing cell lineages; on the contrary, if ${S}_{\mathrm{K}\mathrm{L}}^{(1)}[D]<{S}_{\mathrm{K}\mathrm{L}}^{(2)}[D]$, the population harbors small populations of fastgrowing cell lineages. Untangling the linkage between the structures of growth heterogeneity and their adaptability would help us understand the adaptive strategies of various organisms.
In general, heredity is also crucial for the growth and evolution of a population. The role of the heredity of a particular trait might be unravelled by taking the correlation length as a lineage trait $X$ and quantifying its selection strength. Since the modes of heredity can also be important targets of natural selection (Rivoire and Leibler, 2014), such measurements might provide insights into the evolution of heredity.
We remark that the distribution of interdivision time (generation time) influences the longterm growth rate, as demonstrated by the analytical model in Appendix 2. Therefore, statistical properties of generation time, such as distribution shapes and transgenerational correlations, can contribute to organisms’ evolutionary success by constantly introducing selection within a population. Unlike the central limit theorem, the contributions of higherorder cumulants can remain even in the longterm limit. Importantly, even when cell division processes seem purely stochastic, different states in some traits might underlie these variations in generation times. In such cases, $h(x)$ and ${S}_{\mathrm{rel}}[X]$ for these traits can still unravel the correlations between the trait values and fitness.
This framework is applicable even to cell populations growing under nonconstant environmental conditions. We indeed utilized this framework to analyze the regrowth of growtharrested cells from the stationary phase conditions. The selection strength contributions to population growth, ${S}_{\mathrm{KL}}^{(1)}[D]/\tau \mathrm{\Lambda}$, were below 10% in most cases under constant growth conditions. Nevertheless, it became over 70% in the regrowth of E. coli from the late stationary phase. While increased selection in nonconstant environments may not be surprising itself, it is intriguing to ask how its contribution changes quantitatively depending on the conditions of environmental changes, such as nutrient upshift and downshift. The selection strength contribution in the regrowth from the early stationary phase was only 5%. This result clearly shows that how strongly selection acts in regrowing processes depends on stationary phase incubation durations. However, we also remark that the differences in the selection strength values depend on the time window and might be valid only in this time scale. Clarifying the differences in the selection strength in longer time scales requires the detail of their lag time distributions, which we did not measure in this study.
We identified the cellular populations in which selection acts to increase fitness variance in the retrospective statistics compared with the chronological statistics (Figures 5F and 6G and Figure 5—figure supplement 2). When a decrease in fitness variance by selection is mentioned in evolutionary biology, an upper bound and inheritance of fitness across the generations of individuals are usually assumed. In such circumstances, selection drives the fitness distribution toward the maximum value, and the selection eventually causes fitness variance to decrease. However, even in this process, a decrease is not assured for every step; whether selection reduces fitness variance at each step depends on the fitness distribution at that time. Likewise, whether the fitness variance increases or decreases in the retrospective distribution depends on the shape of the fitness distribution before selection, that is, chronological distribution. Such conditions are graphically recognized by the downward convexity of ${K}_{D}^{\prime}(\xi )$ (Figure 3). When the fourth or higher order fitness cumulants are negligible, the convexity of ${K}_{D}^{\prime}(\xi )$ is determined primarily by the skewness of ${Q}_{\mathrm{cl}}(d)$; positive skew of ${Q}_{\mathrm{cl}}(d)$ with a long right tail makes ${K}_{D}^{\prime}(\xi )$ convex downward and $\mathrm{Var}{[\stackrel{~}{h}(D)]}_{\mathrm{rs}}$ greater than $\mathrm{Var}{[\stackrel{~}{h}(D)]}_{\mathrm{cl}}$. This consequence is intuitively understandable since the right tail of ${Q}_{\mathrm{cl}}(d)$ is accentuated in proportion to ${e}^{D}$ by selection, which leads to greater variance of ${Q}_{\mathrm{rs}}(d)$. On the other hand, when the skew is negative with the long left tail, the effect of applying ${e}^{D}$ is to diminish the tail and compress the distribution toward the fittest lineages. It is of note that greater fitness variance in the retrospective statistics is possible even in the longterm limit, as demonstrated by the model in Appendix 2.
We showed that division count heterogeneity among cellular lineages has dual facets: increasing population growth rate while sensitizing populations to perturbations. These two effects are quantitatively represented by ${S}_{\mathrm{KL}}^{(1)}[D]/\tau \mathrm{\Lambda}$ and ${S}_{\mathrm{KL}}^{(2)}[D]/\tau \mathrm{\Lambda}$, respectively. Therefore, the difference between these selection strength measures gauges which aspect of growth heterogeneity is more significant in the population. Even though ${S}_{\mathrm{KL}}^{(1)}[X]$ and ${S}_{\mathrm{KL}}^{(2)}[X]$ are different in general, the analysis revealed that they were nearly identical in most of the cellular populations growing at constant rates (Figure 5). This result might suggest that, from a practical viewpoint, the contribution of higherorder cumulants becomes negligible under steady growth conditions, and the significant difference between ${S}_{\mathrm{K}\mathrm{L}}^{(1)}[X]$ and ${S}_{\mathrm{KL}}^{(2)}[X]$ could be used as a probe for the nonstationarity of the population growth. This speculation must be examined experimentally using various organisms and cell types across diverse environmental conditions.
This framework is premised on complete lineage tree information. However, many methods of singlecell measurements continuously exclude cells from observation areas and provide only a part of the tree information. Therefore, extending this framework so that one can infer both chronological and retrospective probabilities from incomplete tree information is an essential future research direction. In this study, we calculated the fitness landscapes and selection strength measures for the cell lineage data obtained with the mother machine devices, assuming that these cell lineages would follow the chronological statistics. Such a simple approach is not yet available for larger scale lineage tree data obtainable with the other singlecell measurement devices such as dynamics cytometer (Hashimoto et al., 2016) and chemoflux (Lambert et al., 2014). Furthermore, it has been shown that the inference precision of population growth rate has nonmonotonic dependence on the length of cell lineages obtained with mother machine devices (Levien et al., 2020). Even though the difficulties to overcome are present, a comprehensive framework may permit a unified treatment of cellular lineage data obtained using various singlecell measurement methods.
Phenotypic heterogeneity is widely observed in diverse cellular systems, including both prokaryotic and eukaryotic cells. It is often considered that phenotypic heterogeneity allows bethedging against unpredictable environments and promotes the survival of cellular population (Kussell and Leibler, 2005). However, quantitative evaluation of correlations between the traits of interest and fitness is usually an intricate problem. The cell lineage statistical framework described in this study offers a straightforward procedure applicable to any cellular genealogical data, which are now becoming increasingly available for various biological phenomena, including cancer metastasis (Quinn et al., 2021) and stem cell differentiation (Filipczyk et al., 2015; Frieda et al., 2017; Chow et al., 2021). Another important advantage of this framework is that it allows decomposing a population growth rate into chronological fitness and selection strength. It is thus intriguing to apply this framework to longterm evolutionary dynamics and quantify how the contributions of chronological mean fitness and selection underlie the transitions of population growth rate. Such analysis might clarify the crucial roles of phenotypic heterogeneity in facilitating evolution.
Materials and methods
Microfabrication of microchamber array
We constructed and used a microchamber array for conducting singlecell timelapse 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 exponentialphase 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 stationaryphaseregrowth 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 (DDB201TW, Neoark) on mask blanks (CBL4006DuAZP, CLEAN SURFACE TECHNOLOGY). The photoresist on mask blanks was developed in NMD3 (Tokyo Ohka Kogyo). The uncovered chromium (Cr)layer was removed in MPME30 (DNP Fine Chemicals), and the remaining photoresist was removed by acetone. Lastly, the slide was rinsed in MilliQ water and airdried.
The microchamber array was created in glass coverslips by chemical etching. First, we coated a 1,000angstrom Crlayer on a clean coverslip (NEO Micro glass, No. 1., 24 mm × 60 mm, Matsunami) by evaporative deposition and AZP1350 (AZ Electronic Materials) by spincoating on the Crlayer. We transferred the photomask patterns using a mask aligner (MA20, Mikasa). After developing the photoresist in NMD3 and the Crlayer in MPME30, the coverslip was soaked in buffered hydrofluoric acid solution (110BHF, 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 Crlayer were removed by acetone and MPME30, respectively.
Fabrication of PDMS pad
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 bubbletrap groove, which prevents interference with brightfield microscopic imaging by air bubbles in flowing media.
To create a mold for the bubbletrap groove, we spincoated SU8 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 SU8 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 SU8 photoresist was developed in the SU8 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 SU8 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 10cm silicone tubes (SR1554, 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 singlecell measurements.
Chemical decoration of coverslip and cellulose membrane
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(2aminoethylaminopropyl)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 NHSLCLCBiotin (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 streptavidindecorated 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 Pretreated 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_{4} 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.
E. coli strains
We used two E. coli strains: MG1655 and MG1655 F3 rpoSmcherry (MG1655 ΔfliCΔfimAΔflu rpoSmcherry/pUA66PrplSgfp). MG1655 was used in the regrowth experiment from the stationary phases (Figure 6). MG1655 F3 rpoSmcherry was used for analyzing the growth in constant environments (Figures 5 and 7). In MG1655 F3 rpoSmcherry, the three genes, fliC, fimA, and flu, were deleted, and mcherry gene was inserted downstream of rpoS gene to express RpoSmCherry translational fusion protein. This strain also expresses green fluorescent protein (GFP) from a lowcopy plasmid, pUA66PrplSgfp, taken from a comprehensive library of fluorescent transcriptional reporters (Zaslaver et al., 2006).
Culture conditions and sample preparation (exponential growth)
We used MG1655 F3 rpoSmcherry 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 singlecell observation, we first inoculated a glycerol stock into a 3ml culture medium and incubated it with shaking overnight under the same conditions of culture medium and temperature as those used in the timelapse measurement. 30 μl of the overnight culture was inoculated in a 3ml fresh medium and incubated with shaking until the optical density at $\lambda $ = 600 nm reaches 0.10.3. This exponentialphase culture was diluted to OD_{600} = 0.05, and 0.5 μl of the diluted cell suspension was spotted on the microchamber array on a biotindecorated coverslip. A 5mm × 5mm streptavidindecorated 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 streptavidinbiotin bonding was formed between the coverslip and the membrane. After 5min incubation, the agar pad was removed, and the PDMS pad for medium perfusion was attached on the coverslip via a squareframe twosided seal (FrameSeal Incubation Chambers, Biorad). We immediately filled the device with the fresh medium and connected it to a syringe pump on the microscope stage.
Culture conditions and sample preparation (regrowth from stationary phases)
We used E. coli MG1655 strain and cultured the cells in LuriaBertani (LB) medium at 37°C. To prepare the cells for the timelapse 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_{600} = 0.005 and again cultured with shaking as a preculture. For preparing the earlystationaryphase conditioned medium, 7 ml preculture cell suspension at 8 hr (OD_{600} ≈ 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 timelapse observation, a 10μl preculture cell suspension at 8 hr was mixed with 240 μl earlystationaryphase conditioned medium. One μl of this diluted cell suspension was placed on the microchamber array on a biotindecorated glass coverslip. A 5mm × 5mm streptavidindecorated 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 5min incubation, the conditioned medium agar pad was removed, and the PDMS pad for medium perfusion was attached on the coverslip via a squareframe twosided 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 timelapse observation. The conditioned medium was flushed away immediately before starting the timelapse 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 timelapse 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 24hr preculture cell suspension (OD_{600} ≈ 3.0).
Timelapse measurements and image analysis
We used Nikon TiE inverted microscope equipped with Plan Apo $\lambda $ 100× phase contrast objective (NA1.45), ORCAR2 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 2530 microchambers in parallel in one measurement and acquired the phasecontrast, RpoSmCherry fluorescence, and GFP fluorescence images from each position with a 3min interval. We repeated the timelapse measurement for each culture condition three times. In the regrowth experiment from the stationary phases, we monitored 150250 microchambers in parallel with a 3min interval and acquired only phasecontrast images.
We analyzed the timelapse images by ImageJ (Schneider et al., 2012). We extracted the information of cell size (projected cell area), RpoSmCherry 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.
Data analysis
Distributions and selection strength measures for division count
We calculated the distributions and selection strength measures of $D$ as follows. With the list of division counts $\{D\}$ for each lineage $\sigma $, the chronological and retrospective probabilities were evaluated as ${P}_{\mathrm{cl}}(\sigma )={2}^{D(\sigma )}/{N}_{0}$ and ${P}_{\mathrm{rs}}(\sigma )=1/{N}_{\tau}$, respectively, where N_{0} is the number of cells at $t=0$ and ${N}_{\tau}$ is that at $t=\tau $. From these probabilities, the chronological and retrospective distributions of $D$ were obtained by summing the lineage probabilities for each division count, that is,
The selection strength measures, ${S}_{\mathrm{KL}}^{(1)}[D]$ and ${S}_{\mathrm{KL}}^{(2)}[D]$, were calculated as
where ${D}_{\mathrm{supp}}$ is the support of both chronological and retrospective probabilities with respect to $D$, which is common between the two probabilities.
Distributions and selection strength measures for timeaveraged fluorescence intensity of RpoSmCherry and GFP
We obtained the mean fluorescence intensity of RpoSmCherry and GFP along with the genealogical trees in the timelapse measurements of E. coli MG1655 F3 rpoSmcherry strain. We analyzed the timeaveraged fluorescence intensity of RpoSmCherry and GFP as a lineage trait $X$ and evaluated their distributions, fitness landscapes, and selection strength measures (Figure 7). For each cell lineage, the timeaveraged fluorescence intensity was calculated as
where ${t}_{i}={t}_{\mathrm{start}}+i\mathrm{\Delta}t$ min (t_{start} is the start time of the cell lineage; $\mathrm{\Delta}t=3$ min is the timelapse interval), and ${x}_{\sigma}({t}_{i})$ is the mean fluorescence intensity at time t_{i}.
Generally, bin sizes for the fluorescence intensity affect the selection strength values. However, one can usually find the ranges of bin sizes where the results are relatively insensitive to the choice (Nozoe et al., 2017). Following an empirical rule, we set the bin width $\mathrm{\Delta}X$ to
where $\mathrm{I}\mathrm{Q}\mathrm{R}(X)$ is the interquartile range of the set of $X(\sigma )$ from all the cell lineages. Then, the interval was defined as ${I}_{x,\Delta X}=\left[x\frac{\Delta X}{2},x+\frac{\Delta X}{2}\right]$ for $x=\mathrm{min}(\{X\}),\mathrm{min}(\{X\})+\mathrm{\Delta}X,\mathrm{\cdots},\mathrm{min}(\{X\})+(L1)\mathrm{\Delta}X$, where $L$ is the number of total bins given by $L=\lfloor \frac{\mathrm{max}(\{X\})\mathrm{min}(\{X\})}{\mathrm{\Delta}X}\rfloor +2$.
We calculated the chronological and retrospective probability distributions of $X$ by
$h(x)$ The fitness landscape was evaluated by
The selection strength measures were evaluated by
Cumulant generating functions and cumulants
To plot the differential of the cumulant generating functions in Figure 5FH, we evaluated $K}_{D}^{\mathrm{\prime}}(\xi )=\frac{\sum _{d\in {D}_{\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}}}(d\mathrm{ln}2){2}^{\xi d}{Q}_{\mathrm{c}\mathrm{l}}(d)}{\sum _{d\in {D}_{\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}}}{2}^{\xi d}{Q}_{\mathrm{c}\mathrm{l}}(d)$ by changing $\xi $ from 0 to 1 with the step size 0.01.
We calculated the cumulative contributions of fitness cumulants to the population growth ${W}_{n}^{(X)}$ (Figures 5A—7FH) using a julia package, JuliaDiff/TaylorSeries.jl (Benet and Sanders, 2019; Benet and Sanders, 2021).
Error estimations by resampling method
To evaluate the error ranges of the quantities calculated in the analysis, we created 20,000 randomly resampled datasets for each condition and reported the means and two standard deviation ranges in the results.
For the datasets of colony growth (E. coli and M. smegmatis), ${N}_{\tau}$ lineages were randomly sampled with replacement according to the probability weight ${P}_{\mathrm{rs}}(\sigma )$ for each resampled dataset. In each resampled dataset, the initial number of cells was estimated as ${\widehat{N}}_{0}={\sum}_{\sigma \in {\{\sigma \}}_{\mathrm{resampled}}}{2}^{D(\sigma )}$.
For the datasets taken using the mother machines (S. pombe and L1210), we randomly sampled N_{0} lineages with an equal weight, which corresponds to the chronological probability in this setting. ${N}_{\tau}$ was estimated as ${\widehat{N}}_{\tau}={\sum}_{\sigma \in {\{\sigma \}}_{\mathrm{resampled}}}{2}^{D(\sigma )}$.
Simulating the effect of cell removal on population growth rates
We simulated cell population growth with cell removal using a custom C script. The gamma distributions were adopted as generation time distributions. We assigned the shape parameter to $k=$ 1, 2, or 5 and the scale parameter to $\theta ={2}^{1/k}1$. The perturbation strength $\u03f5$ was changed from 0 to 0.2 with the interval 0.01.
As a prerun, we started a simulation from a newborn cell and assigned its generation time randomly according to a predefined gamma probability distribution. We assumed that this cell divided into two daughter cells at the end of the generation. Each daughter cell was removed with probability $1{2}^{\u03f5}$ and assigned with generation time from the same predefined probability distribution if it escaped removal. Repeating this procedure, we let the population grow until all of the remaining cell lineages in the population exceed the maximum duration ${T}_{\mathrm{max}}=8.0$. The time to the next division of each cell lineage at ${T}_{\mathrm{max}}$ was exported as the first division time in the main simulation. This prerun was repeated 1000 cycles to export a sufficiently sizable list of first division times.
In the main simulation, we started from a progenitor cell with its division time randomly assigned from the first division time list exported in the prerum. For the daughter cells born from the first divisions and their descendants, the assignment of generation time and the cell removal were done as in the prerun. We stopped further production of daughter cells in each lineage if it exceeded ${T}_{\mathrm{max}}=8.0$. We repeated this main simulation 1,000 cycles starting from different progenitor cells. The number of cell divisions in each cell lineage until ${T}_{\mathrm{max}}$ was exported for analysis.
We calculated the population growth rate at each perturbation strength as
where $N({T}_{\mathrm{max}},\u03f5)$ is the number of cell lineages at ${T}_{\mathrm{max}}$ when the perturbation strength was $\u03f5$. The chronological and retrospective mean fitness of division count without cell removal was calculated as
When simulating the cell population with motherdaughter correlation time, we randomly assigned the generation time from the gamma probability distribution with its shape parameter $\frac{r{\tau}_{m}/\theta +k(1r)}{1{r}^{2}}$ and scale parameter $(1{r}^{2})\theta $, where ${\tau}_{m}$ is the generation time of the mother cell, $r$ is the correlation coefficient of generation time between neighboring generations. The stationary distribution of this transition probability approximates the gamma distribution with shape parameter $k$ and scale parameter $\theta $ to good precision with identical first and secondorder moments irrespective of the parameters $k$, $\theta $, and $r$. In Figure 4—figure supplement 1, we fixed $k=2$ and $\theta =\sqrt{2}1$ and set $r$ to 0, 0.2, 0.4, or 0.6.
Data and code availability
The raw data obtained in this study, the Matlab codes for data analysis, and the C code for simulation have been deposited in Github repositories (https://github.com/Wakamotolab/LineageAnalysis, (copy archived at swh:1:rev:1865d167f1c24625c98d3c493a9a180b1aa2035d; Yamauchi, 2021), https://github.com/Wakamotolab/LineageAnalysisJulia, (copy archived at swh:1:rev:e22fbce8a713582a18fbe2bcc57dc9078090f121; Nozoe and Wakamoto, 2021) and https://github.com/Wakamotolab/LineageSimulation, (copy archived at swh:1:rev:ef1166620396835168ca9061851898993a091976; Wakamoto, 2021).
Appendix 1
Analytical calculations of fitness measures, selection strength, and the cumulants of a fitness landscape
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; $\varphi <1$ for binomial; and $\varphi >1$ for negative binomial (Appendix 1—figure 1A). Then, ${K}_{D}(\xi )$ for these distributions have a closed form
(Appendix 3). We then immediately obtain
Since ${lim}_{\varphi \to 2}{K}_{D}\left(1\right)=\mathrm{\infty}$, $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 $W}_{n}^{(D)$ (Equation 14). Plotting ${W}_{n}^{(D)}$ shows that the contribution of higherorder 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 $\varphi =0,{\varphi}_{0}(=0.5857)$; ${S}_{\mathrm{K}\mathrm{L}}^{\left(1\right)}\left[D\right]>{S}_{\mathrm{K}\mathrm{L}}^{\left(2\right)}\left[D\right]$ if $0<\varphi <{\varphi}_{0}$ and ${S}_{\mathrm{K}\mathrm{L}}^{\left(1\right)}\left[D\right]<{S}_{\mathrm{K}\mathrm{L}}^{\left(2\right)}\left[D\right]$ if ${\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)$.
Appendix 2
Longterm limit for gammadistributed uncorrelated generation times
To understand how inherent stochasticity affect longterm 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 socalled EulerLotka 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.,
where $\alpha >0$ is a shape parameter; and $\theta >0$ is a scale parameter. In this case, the Malthusian parameter is
The probability distribution of division count ${Q}_{\mathrm{cl}}(D)$, in this case, is known as gamma count distribution (Winkelmann, 1995). Though any closedform expression of the corresponding cumulant generating function is not known, it has a simple limiting form for $\tau \to \mathrm{\infty}$ as shown below. We define the rescaled cumulant generating function by
Since ${\stackrel{~}{K}}_{D}\left(\xi \right)$ represents the population growth rate, or Malthusian parameter with the mean number of offspring $z={2}^{\xi}$, we have
When $g$ is a gamma distribution with a shape parameter $\alpha $ and a scale parameter $\theta $, we obtain
and
Note that $\alpha =1$ corresponds to the case where division counts follow the Poisson distribution with mean ${\theta}^{1}$. The scaled key quantities derived from ${\stackrel{~}{K}}_{D}(\xi )$ are as follows.
and
Hence,
and
${\stackrel{~}{S}}_{\mathrm{K}\mathrm{L}}^{\left(2\right)}\left[D\right]>{\stackrel{~}{S}}_{\mathrm{K}\mathrm{L}}^{\left(1\right)}\left[D\right]$ is always true for $0<\alpha <\mathrm{\infty}$ because
where $\gamma ={\alpha}^{1}\mathrm{ln}2$ and the inequality ${e}^{\gamma}>1+\gamma$ ($\gamma >0$) are used.
Since the Taylor expansion of ${\stackrel{~}{K}}_{D}\left(\xi \right)$ at $\xi =0$ is
the timescaled $n$th order fitness cumulant is
Therefore,
These results show that, unlike the central limit theorem, higherorder cumulants remain even in the longterm limit. Selection strength also remains in the longterm limit, which means that inherent stochasticity of generation times continuously introduces selection within a cellular population. Importantly, the timescaled cumulants and the selection strength depend on $\alpha $. Therefore, the shape of generation time distributions influences the longterm population growth rate and selection. Since ${\stackrel{~}{S}}_{\mathrm{KL}}^{\left(2\right)}\left[D\right]$ is always greater than ${\stackrel{~}{S}}_{\mathrm{KL}}^{\left(1\right)}\left[D\right]$, the fitness variance is larger in the retrospective distribution than in the chronological distribution.
Appendix 3
The properties of the selection strength of division count
Below we derive several important properties of the selection strength of division count. We focus on the selection strength measure ${S}_{\mathrm{KL}}^{(1)}$ and write it as $S$ this section for conciseness. However, the conclusions are likewise valid for ${S}_{\mathrm{JF}}$ and ${S}_{\mathrm{KL}}^{(2)}$.
The most detailed description of cellular lineage statistics is based on individual lineages $\sigma $. From the definitions of ${P}_{\mathrm{cl}}(\sigma )$ and ${P}_{\mathrm{rs}}(\sigma )$ in the main text, the relation
is held. We define the selection strength of cellular lineages as
where ${\u27e8D(\sigma )\mathrm{ln}2\u27e9}_{\mathrm{cl}}={\sum}_{\sigma}(D(\sigma )\mathrm{ln}2){P}_{\mathrm{cl}}(\sigma )$
From the definition of fitness landscape (Equation 1),
On the other hand,
This proves that the chronological mean fitness of cellular lineages equals the chronological mean fitness of division count.
Since $S[D]=\tau \mathrm{\Lambda}{\u27e8\stackrel{~}{h}(D)\u27e9}_{\mathrm{cl}}$ and $S[\sigma ]=\tau \mathrm{\Lambda}{\u27e8D(\sigma )\mathrm{ln}2\u27e9}_{\mathrm{cl}}$ (Equations 3; 58),
is also held. This result shows that the selection strength of $D$ is equivalent to the selection strength of cellular lineages despite $D$ being a coarsegrained lineage trait.
Another important property of $S[D]$ is that it sets the maximum bound for the selection strength of any lineage traits. Now we consider the joint probability distributions of $D$ and lineage trait $X$, which we write ${Q}_{\mathrm{cl}}(d,x)$ and ${Q}_{\mathrm{rs}}(d,x)$. We define the joint selection strength as
Using ${Q}_{\mathrm{cl}}(d,x)={Q}_{\mathrm{cl}}(dx){Q}_{\mathrm{cl}}(x)$ and ${Q}_{\mathrm{rs}}(d,x)={Q}_{\mathrm{rs}}(dx){Q}_{\mathrm{rs}}(x)$,
where $S[DX]:=\sum _{d}\sum _{x}{Q}_{\mathrm{c}\mathrm{l}}(d,x)\mathrm{ln}\frac{{Q}_{\mathrm{c}\mathrm{l}}(dx)}{{Q}_{\mathrm{r}\mathrm{s}}(dx)}$, and we used $\sum _{d}{Q}_{\mathrm{c}\mathrm{l}}(dx)=1$.
Likewise, $S[D,X]$ can also be decomposed as
However, $S[XD]=0$ because
and
from (Equation 1) and (Equation 65). This leads to
from (Equation 63) and (Equation 64). Furthermore, $S[DX]\ge 0$ from Jensen’s inequality. Thus,
The equality is held when $D$ is a deterministic function of $X$. This inequality shows that $S[D]$ ($=S[\sigma ]$) sets the maximum bound for the selection strength of any lineage trait $X$.
The cumulant generating function ${K}_{X}(\xi )$ provides both chronological and retrospective fitness cumulants
In the main text, we introduced the cumulant generating function of $h(x)$ with respect to the chronological distribution ${Q}_{\mathrm{cl}}(x)$,
This function can also be written as
when the fitness cumulants ${\kappa}_{n}^{(X)}$ are all finite, and the Taylor expansion converges at $\xi $. Also,
Below we prove that ${K}_{X}(\xi )$ also gives the fitness cumulants on the retrospective distributions.
We define a cumulant generating function on the retrospective probability as
This function can be expanded by the fitness cumulants of the retrospective statistics ${\rho}_{n}^{(X)}$ as
Therefore,
For example, ${\rho}_{1}^{(X)}={\u27e8h(X)\u27e9}_{\mathrm{rs}}$ and ${\rho}_{2}^{(X)}=\mathrm{Var}{[h(X)]}_{\mathrm{rs}}={\u27e8h{(X)}^{2}\u27e9}_{\mathrm{rs}}{\u27e8h(X)\u27e9}_{\mathrm{rs}}^{2}$.
Inserting ${Q}_{\mathrm{rs}}(x)={e}^{h(x)\tau \mathrm{\Lambda}}{Q}_{\mathrm{cl}}(x)$ into (Equation 72),
Hence,
for $n\ge 1$. This relation proves that evaluating $\frac{{d}^{n}{K}_{X}(\xi )}{d{\xi}^{n}}$ at $\xi =1$ gives the $n$th order fitness cumulant on the retrospective statistics; i.e.,
Furthermore, this leads to
from (Equation 70) and (Equation 77). Similarly, evaluating (Equation 76) at $\xi =1$ gives
Analogously to (Equation 12), we can also expand the population growth rate in terms of the retrospective cumulants, by evaluating (Equation 75) at $\xi =1$,
For example, when the fitness distribution is Gaussian for the chronological statistics,
since ${\kappa}_{n}^{X}=0$ for $\forall n\ge 3$.
These results confirm that the function ${K}_{X}(\xi )$ contains the information of both chronological and retrospective statistics.
Relationships between fitness cumulants and selection strength measures
In the main text, we have shown that the selection strength ${S}_{\mathrm{KL}}^{(1)}[X]$ corresponds to the contribution of the second or higherorder fitness cumulants to population growth, i.e.,
or alternatively, by substituting (Equation 79) and (Equation 80) we obtain
Similar expressions can also be found for ${S}_{\mathrm{KL}}^{(2)}[X]$. Since ${S}_{\mathrm{KL}}^{(2)}[X]={\u27e8h(X)\u27e9}_{\mathrm{rs}}\tau \mathrm{\Lambda}$ (Equation 4), substituting (Equation 80) yields
or alternatively, by substituting (Equation 78) and (Equation 12) we obtain
These show that both of ${S}_{\mathrm{KL}}^{(1)}[X]$ and ${S}_{\mathrm{KL}}^{(2)}[X]$ can be expanded by the chronological or retrospective fitness cumulants.
The difference between these two selection strength measures is
from (Equation 83) to (Equation 86). Thus, it depends only on the third or higherorder fitness cumulants.
Finally, another selection strength measure ${S}_{\mathrm{JF}}[X]$ can also be expanded by the fitness cumulants as
from (Equation 83) to (Equation 86). When the chronological fitness distribution is Gaussian (${\kappa}_{n}^{(X)}=0$ for $\forall n\ge 3$),
Analytical calculations of ${K}_{D}(\xi )$ and related relations given specific form of division count distributions
Here we derive (Equations 35–38) in the main text. We begin with the case where ${Q}_{\mathrm{cl}}(D)$ follows a Poisson distribution. Let $\overline{D}$ denote the chronological mean division count.
By the definition of ${K}_{D}(\xi )$,
By the Taylor expansion of ${2}^{\xi}={e}^{\xi \mathrm{ln}2}$, the $n$th order cumulant is ${\kappa}_{n}^{(D)}=\overline{D}{\left(\mathrm{ln}2\right)}^{n}$. Since
we derive
For example, ${W}_{1}^{(D)}=0.693$, ${W}_{2}^{(D)}=0.933$, ${W}_{3}^{(D)}=0.988$, and ${W}_{4}^{(D)}=0.998$. The first order derivative of ${K}_{D}(\xi )$ is
and thereby we have
and
Next we derive ${K}_{D}(\xi )$ for binomial and negative binomial distributions. Let $\overline{D}$ and $\overline{D}\varphi $ denote the mean and the variance of ${Q}_{\text{cl}}\left(D\right)$. When ${Q}_{\text{cl}}\left(D\right)$ is binomial,
where ${D}_{\mathrm{max}}$ and $p$ satisfy $\overline{D}={D}_{\mathrm{max}}p$ and $\overline{D}\varphi ={D}_{\mathrm{max}}p(1p)$; namely ${D}_{\mathrm{max}}=\overline{D}/\left(1\varphi \right)$ and $p=1\varphi $. Therefore,
When ${Q}_{\text{cl}}\left(D\right)$ is negative binomial,
where $\alpha $ and $p$ satisfy $\overline{D}=\alpha p/(1p)$ and $\overline{D}\varphi =\alpha p/{(1p)}^{2}$; namely $\alpha =\overline{D}/\left(\varphi 1\right)$ and $p=1{\varphi}^{1}$. Therefore,
(Equation 102) is exactly the same as (Equation 100) as the function of $\overline{D},\varphi ,$ and $\xi $. In addition, (Equation 91) is the limiting form of (Equation 100) and (Equation 102) as $\varphi \to 1$. Thus, (Equation 35) in the main text represents ${K}_{D}(\xi )$ for Poisson, binomial or negative binomial ${Q}_{\mathrm{cl}}(D)$.
The Taylor expansion of (Equation 35) is obtained as follows:
where
For the first five terms, for example, we have
${\kappa}_{n}^{(D)}=\overline{D}{c}_{n}(\varphi ){(\mathrm{ln}2)}^{n}$ gives the $n$th order cumulant.
The first order derivative of (Equation 35) is
and thereby we obtain
and
(Equation 109) and (Equation 110) equal if and only if
This equation has two roots $\varphi =0$ and $\varphi ={\varphi}_{0}=0.5857\mathrm{\dots}$ and LHS >RHS if and only if $0<\varphi <{\varphi}_{0}$.
Data availability
All data generated or analyzed during this study and the Matlab codes for data analysis have been deposited in a GitHub repository (https://github.com/Wakamotolab/LineageAnalysis; copy archived at swh:1:rev:1865d167f1c24625c98d3c493a9a180b1aa2035d).

Dryad Digital RepositoryData from: Inferring fitness landscapes and selection on phenotypic states from singlecell genealogical data.https://doi.org/10.5061/dryad.4539d

Dryad Digital RepositoryData from: Aging, mortality, and the fast growth tradeoff of Schizosaccharomyces pombe.https://doi.org/10.5061/dryad.s2t5t

Dryad Digital RepositoryData from: Intrinsic growth heterogeneity of mouse leukemia cells underlies differential susceptibility to a growthinhibiting anticancer drug.https://doi.org/10.5061/dryad.80gb5mkpr
References

Bacterial persistence as a phenotypic switchScience 305:1622–1625.https://doi.org/10.1126/science.1099390

The rposmediated general stress response in Escherichia coliAnnual Review of Microbiology 65:189–213.https://doi.org/10.1146/annurevmicro090110102946

TaylorSeries.jl: taylor expansions in one and several variables in juliaJournal of Open Source Software 4:1043.https://doi.org/10.21105/joss.01043

Advanced methods of microscope control using μmanager softwareJournal of Biological Methods 1:e10.https://doi.org/10.14440/jbm.2014.36

Stochastic gene expression in a single cellScience 297:1183–1186.https://doi.org/10.1126/science.1070919

Network plasticity of pluripotency transcription factors in embryonic stem cellsNature Cell Biology 17:1235–1246.https://doi.org/10.1038/ncb3237

BookThe Genetical Theory of Natural SelectionOxford: Oxford University Press.https://doi.org/10.5962/bhl.title.27468

Natural selection. V. how to read the fundamental equations of evolutionary change in terms of information theoryJournal of Evolutionary Biology 25:2377–2396.https://doi.org/10.1111/jeb.12010

Universal constraints on selection strength in lineage treesPhysical Review Research 3:023187.https://doi.org/10.1103/PhysRevResearch.3.023187

Fundamental principles in bacterial physiologyhistory, recent progress, and the future with focus on cell size control: a reviewReports on Progress in Physics. Physical Society 81:056601.https://doi.org/10.1088/13616633/aaa628

The growth rate of individual bacterial cellsJournal of Bacteriology 23:147–153.https://doi.org/10.1128/jb.23.2.147153.1932

Fluctuation relations of fitness and information in population dynamicsPhysical Review Letters 115:238102.https://doi.org/10.1103/PhysRevLett.115.238102

A theory for the age and generation time distribution of a microbial populationJournal of Mathematical Biology 1:17–36.https://doi.org/10.1007/BF02339486

Large deviation principle linking lineage statistics to fitness in microbial populationsPhysical Review Letters 125:048102.https://doi.org/10.1103/PhysRevLett.125.048102

Cell cycle heritability and localization phase transition in growing populationsPhysical Review Letters 125:268103.https://doi.org/10.1103/PhysRevLett.125.268103

SoftwareLineageAnalysisjulia, version swh:1:rev:e22fbce8a713582a18fbe2bcc57dc9078090f121Software Heritage.

Growth rate and generation time of bacteria, with special reference to continuous cultureJournal of General Microbiology 15:492–511.https://doi.org/10.1099/00221287153492

Ergodicity, hidden bias and the growth rate gainPhysical Biology 15:036006.https://doi.org/10.1088/14783975/aab0e6

Nih image to imagej: 25 years of image analysisNature Methods 9:671–675.https://doi.org/10.1038/nmeth.2089

Optimal lineage principle for agestructured populationsEvolution; International Journal of Organic Evolution 66:115–134.https://doi.org/10.1111/j.15585646.2011.01418.x

SoftwareLineageSimulation, version swh:1:rev:ef1166620396835168ca9061851898993a091976Software Heritage.

Robust growth of Escherichia coliCurrent Biology 20:1099–1103.https://doi.org/10.1016/j.cub.2010.04.045

Duration dependence and dispersion in countdata modelsJournal of Business & Economic Statistics 13:467–474.https://doi.org/10.1080/07350015.1995.10524620

SoftwareLineageAnalysis, version swh:1:rev:1865d167f1c24625c98d3c493a9a180b1aa2035dSoftware Heritage.
Article and author information
Author details
Funding
Japan Science and Technology Agency (JPMJCR1927)
 Yuichi Wakamoto
Japan Science and Technology Agency (JPMJER1902)
 Yuichi Wakamoto
National Institute of General Medical Sciences (R01GM097356)
 Edo Kussell
Japan Society for the Promotion of Science (17H06389)
 Yuichi Wakamoto
Japan Society for the Promotion of Science (19H03216)
 Yuichi Wakamoto
Japan Society for the Promotion of Science (21K20672)
 Takashi Nozoe
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Tetsuya J Kobayashi and the members of the Wakamoto Lab for discussion. This work was supported by JST CREST Grant Number JPMJCR1927 (YW); JST ERATO Grant Number JPMJER1902 (YW); NIH Grant Number R01GM097356 (EK); and Japan Society for the Promotion of Science KAKENHI Grant Number 17H06389 (YW), 19H03216 (YW), and 21K20672 (TN).
Version history
 Preprint posted: July 19, 2021 (view preprint)
 Received: July 26, 2021
 Accepted: October 28, 2022
 Version of Record published: December 6, 2022 (version 1)
Copyright
© 2022, Yamauchi et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 580
 views

 82
 downloads

 3
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
Transcriptomic profiling became a standard approach to quantify a cell state, which led to accumulation of huge amount of public gene expression datasets. However, both reuse of these datasets or analysis of newly generated ones requires significant technical expertise. Here we present Phantasus  a userfriendly webapplication for interactive gene expression analysis which provides a streamlined access to more than 96000 public gene expression datasets, as well as allows analysis of useruploaded datasets. Phantasus integrates an intuitive and highly interactive JavaScriptbased heatmap interface with an ability to run sophisticated Rbased analysis methods. Overall Phantasus allows users to go all the way from loading, normalizing and filtering data to doing differential gene expression and downstream analysis. Phantasus can be accessed online at https://alserglab.wustl.edu/phantasus or can be installed locally from Bioconductor (https://bioconductor.org/packages/phantasus). Phantasus source code is available at https://github.com/ctlab/phantasus under MIT license.

 Computational and Systems Biology
 Evolutionary Biology
A comprehensive census of McrBC systems, among the most common forms of prokaryotic Type IV restriction systems, followed by phylogenetic analysis, reveals their enormous abundance in diverse prokaryotes and a plethora of genomic associations. We focus on a previously uncharacterized branch, which we denote coiledcoil nuclease tandems (CoCoNuTs) for their salient features: the presence of extensive coiledcoil structures and tandem nucleases. The CoCoNuTs alone show extraordinary variety, with three distinct types and multiple subtypes. All CoCoNuTs contain domains predicted to interact with translation system components, such as OBfolds resembling the SmpB protein that binds bacterial transfermessenger RNA (tmRNA), YTHlike domains that might recognize methylated tmRNA, tRNA, or rRNA, and RNAbinding Hsp70 chaperone homologs, along with RNases, such as HEPN domains, all suggesting that the CoCoNuTs target RNA. Many CoCoNuTs might additionally target DNA, via McrC nuclease homologs. Additional restriction systems, such as Type I RM, BREX, and Druantia Type III, are frequently encoded in the same predicted superoperons. In many of these superoperons, CoCoNuTs are likely regulated by cyclic nucleotides, possibly, RNA fragments with cyclic termini, that bind associated CARF (CRISPRAssociated Rossmann Fold) domains. We hypothesize that the CoCoNuTs, together with the ancillary restriction factors, employ an echeloned defense strategy analogous to that of Type III CRISPRCas systems, in which an immune response eliminating virus DNA and/or RNA is launched first, but then, if it fails, an abortive infection response leading to PCD/dormancy via host RNA cleavage takes over.