Ergodicity, Hidden Bias and the Growth Rate Gain
Many single-cell observables are highly heterogeneous. A part of this heterogeneity stems from age-related phenomena: the fact that there is a nonuniform distribution of cells with different ages. This has led to a renewed interest in analytic methodologies including use of the “von Foerster equation” for predicting population growth and cell age distributions. Here we discuss how some of the most popular implementations of this machinery assume a strong condition on the ergodicity of the cell cycle duration ensemble. We show that one common definition for the term ergodicity, “a single individual observed over many generations recapitulates the behavior of the entire ensemble” is implied by the other, “the probability of observing any state is conserved across time and over all individuals” in an ensemble with a fixed number of individuals but that this is not true when the ensemble is growing. We further explore the impact of generational correlations between cell cycle durations on the population growth rate. Finally, we explore the “growth rate gain” - the phenomenon that variations in the cell cycle duration lead to an improved population-level growth rate - in this context. We highlight that, fundamentally, this effect is due to asymmetric division.
In most cell biology experiments, measurements are made on cells during exponential growth, and the results are often very heterogeneous from cell to cell. This has led to the prevalent use of -values and other sophisticated statistical measures to determine if two populations are likely to be different. There are several fundamental reasons for this observed variabilitymetz2014dynamics. ranging from stochastic protein synthesis rates to the natural selection of diverse, adaptable populationsrochman2016grow; bodi2017phenotypic. One of the most important reasons is that cells are at different stages of the cell cycle, and therefore display different cycle-dependent protein expression levels. For cells replicating mitotically, there are more younger cells than older cells present in a random selection of individuals; therefore, the age (as measured from the moment of completion of mitosis) distribution of the cells determines, to a large extent, the statistical distribution of observables. The basic formalism for modeling the age distribution of a population was developed by von Foerster von1960doomsday. Starting from a known cell cycle duration (the time between cell birth and cell division) distribution, the expected age distribution can be explicitly computed powell1956growth; stukalin2013age; however, the von Foerster approach assumes a strong condition for the ergodicity of the cell population. In particular, this formalism implicitly assumes that the cell cycle durations from generation to generation are independently sampled from the same distribution across all individuals. Recent experimental work has shown measurable correlations in cell cycle duration contradicting this assumption for both bacterial wang2010robust; hashimoto2016noise and yeast populations cerulus2016noise. These measurements have been made in a variety of settings including microfluidic devices where bacteria are maintained at constant density and media conditions in channels of around ten to twenty cellswang2010robust and up to one hundred cellshashimoto2016noise, as well as agar sandwiches for yeast where density remains low throughout the experiment, though it is not explicitly controlledcerulus2016noise. Moreover, it has been shown that higher variability in the cell cycle duration distribution leads to higher population-level growth rates cerulus2016noise; hashimoto2016noise, a phenomenon labelled the “growth rate gain”. In this paper, we explicitly examine the relationships between generational correlations within the cell cycle duration distribution, width of the cell cycle duration distribution, and the population-level growth rate. We show that the age distribution obtained from observing the cell cycle duration distribution differs if observations are made of an ensemble of fixed size (e.g. within a microfluidic device where only one of two daughter cells remains in the ensemble under observation after cell division or randomly selected from a dividing population. We discuss how the width of the cell cycle duration distribution affects the generational cell cycle duration correlations predicted using the von Foerster formalism and highlight that, fundamentally, the “growth rate gain” in this context is due to asymmetric division.
Consider the following conservation of population where is the population per unit age at time t, . This states the population of cells at age and time is equal to the population of cells at age and time after removing those that have left the ensemble (due to mitosis, etc.). The loss rate, in units of is notated as . For the purposes of this paper, we may consider “chronological age”, , though in the case where the “age” of interest is cell phase or another non-chronological age, this introduces some added complexity (please seerubinow1968maturity equations 1 and 2). Stated another way, in the case of chronological age, a cell is -minutes older after -minutes of lab-frame time has passed. While this may be straight-forward, it may not always be useful. For example, if you want to compare a fast growing cell with a slow growing cell, both will have the same chronological aging rate ; however, the fast growing cell will progress through the cell cycle much faster. To capture this phenomenon, we may want to examine the cell cycle “phase” rather than the chronological age. In the simpler case, , we find the relatively simple form:
To solve this equation, we can examine the stage when the cell is undergoing exponential growth, where is the population growth rate and is the steady state age distribution. Using this assumed form, we may solve for . To move forward, we need some information about . We are focusing on the case where the cell multiplies through mitosis. In this case, when the death rate is negligible, is just the division rate, and can be obtained directly from the cell cycle duration probability density function, . has been explicitly measured in bacteriawang2010robust; stukalin2013age; rochman2016grow, yeastcerulus2016noise, and mammalian cellsstukalin2013age. We find
is the probability of undergoing mitosis per unit time at age , and is the ratio of the population of cells observed to divide at age over the population of cells which have matured to age without yet dividing. We note that this too is only valid in the chronological case (this is clear from just the units: here , which as defined has units of , takes the units of , ). also appears in the following boundary condition: the number of cells at age zero and time is just twice the number of cells that divided at time : . We may now rewrite this boundary condition in terms of explicitly:
and similarly rewrite the steady state age distribution to yield:
This framework (Eqs. 3 and 4) provides a way to calculate the population growth rate, , and the age distribution, , given only the cell cycle duration distribution, ; however, it is built on some strong assumptions owing to the interpretation of Eq. 1. We have already discussed that this result is only valid for the case of chronological aging, and noted that we assume cell death is negligible. Beyond this, there is a third assumption which leads to some important consequences we want to discuss. We rewrote the division rate in terms of the cell cycle duration distribution, and in doing so, assumed that is the same for every cell, at all times. In other words, consider the transition probability (per unit time) between successive cell cycle durations, where represents the cell cycle duration for the cell of interest during generation . In general, this function may have some dependence on , , and even explicit dependence on time; however, we have strictly assumed:
To see how this assumption impacts the growth rate, let us first take a look at how maps to in exponentially growing ensembles as well as ensembles containing a fixed number of individuals where only one daughter cell remains in the ensemble under observation after division.
Ii Age Distribution Properties
First we will illustrate that, under the condition , the age distribution is often “mean-scalable”. Under many circumstances, belongs to a class of functions which are mean-scalablestukalin2013age (please note thatstukalin2013age does not define the term “mean-scalable” but discusses measured which satisfy the definition below). Let us label the mean of to be , and introduce the variable : we will say is mean-scalable when there exists a scaled function which is conserved across all (and this function satisfies normalization):
Stated another way, suppose we have two functions and with mean values and . Additionally, let us consider the functions = and = . If for all , then the family of functions consisting of and are mean-scalable. We can show that the mean-scalability of confers the same property to . Rewriting Eq. 3 in terms of , , and , the scaled bulk growth rate, we find . Note that this uniquely defines given and that for fixed , increasing decreases . This inverse relationship expresses that cells which take longer to divide result in a slower growing ensemble. Similarly, we may now consider the function where the factor of is introduced again to maintain normalization. Rewriting Eq. 4 yields . We just saw that is completely determined by which implies that is also completely determined by . In other words, whenever is mean-scalable, is also mean-scalable. This property is useful because it has been observed that is often mean-scalable. More specifically, is well represented by a gamma distribution:
where and the coefficient of variation: , across widely differing cell types. In this case and we see that the ensemble is no longer mean-scalable when the (a function of ) changes. We note, as shown in Fig. 1, that even when the for the cell cycle duration distribution does vary and is not mean-scalable, the differences in the resultant scaled age distributions, , are still much smaller than the original cell cycle duration distributions.
As mentioned in the introduction, the von Foerster equation is usually solved after making some important assumptions. One straightforward assumption is that the species studied is undergoing exponential growth ; but as clear as this may seem, it is not valid for some experimentally relevant ensembles. A wide variety of microfluidic devices and extracellular matrix patterns are now available to study bacterial, yeast, and mammalian cell populations at constant density, following a single cell for multiple cycles. In these devices, only one of two daughter cells is maintained in the ensemble after each division. This keeps the total number of cells constant over time and modifies the age distribution. We could modify the von Foerster equation to remove the time dependence and solve using the same method as above; but below, we will utilize a simpler framework. Let us take a moment to look at the age distributions that result from these ensembles.
We begin again with the cell cycle duration distribution . Unlike an exponentially growing ensemble, an ensemble with a fixed number of individuals has very simple lineages: where is the cell cycle duration of generation . Since we are still considering the condition , we may also note that the behavior of each lineage recapitulates that of the entire ensemble (this idea is discussed in detail in the following section on ergodicity). Let us follow generations of a single cell over some window of time and observe just one of the two daughter cells after each division. The probability of observing a cell cycle of duration exactly at an arbitrary time where is simply the duration of the cycle multiplied by the number of times it occurred. . This generalizes to the probability density function:
for observing a cycle of duration . Now we may consider a cycle of duration spanning the window of time . If the cell is observed at an arbitrary time point within that window, the probability the cell will have been within the age range , is simply if and if is smaller. Thus the probability of observing a cell of age at an arbitrary time during a cycle of unknown duration is:
We may note using the same method as above, that when is mean scalable, this ensemble is mean scalable too: . Another simple and straightforward, but useful observation is that the mean age of the fixed ensemble can be calculated with a single integral. Written in the usual way, it is a bit cumbersome to calculate directly: . However, we may more easily write down an expression for the mean age of any lineage, which will be the same as the mean age of the entire ensemble. The mean age of each lineage is simply the average of the mean age of each cycle (the mean age of a cycle of length is simply ):
In the case where is gamma distributed see Eq. 7, this yields a closed form expression:
We see that the mean age remains close to until the gets quite large (since the is rarely greater than 1 for experimentally observed hashimoto2016noise; cerulus2016noise; wang2010robust; rochman2016grow). We may also note that in the delta function limit for the cell cycle duration distribution, is simply uniform. In contrast for the exponentially growing case, in the limit where the tends to zero and is a delta function at , we retrieve . We see in this case, the mean age is:
This is about . In Fig. 2 we show that the fixed population is older than the exponentially growing population when is gamma distributed for all . We also want to emphasize this is true in general, independent of the form of . This can be observed from a comparison of the age distributions: and . The two expressions differ only by the factor which monotonically decreases with age. This means if you examine a snapshot of cells which have been maintained in a population of fixed number (e.g. mother cells in a microfluidic device where only one of two daughter cells remains in the ensemble under observation after cell division), the ensemble will be older than a group of unconstrained cells. Stated another way, if you were to pick two arbitrary cells: one from a population of fixed number and one from an exponentially growing population, the cell from the population of fixed number would probably be the older cell. In Fig. 2, we examine how changes with respect to changing of the cell cycle duration distribution, assuming gamma-distributed , and compare with . We see that there is a stronger dependence on the of for the fixed distributions.
The condition, , which led to the nice properties of the age distributions discussed above is really a statement about the ergodicity of the cell cycle duration distribution which may not be true for the biological system of interest. The following definition is usually used to describe an ergodic system: consider an ensemble of measurements made at time of individual . This ensemble is considered to be ergodic if ; that is, the probability of observing state is the same at anytime and from any individual. To clarify, this work focuses on the case where the state observed is the cell cycle duration. Alternatively, this can be written:
where is the probability of observing state conserved across all individuals at all times. We want to illustrate, for some ensembles, that this condition implies over time the same behavior is recapitulated as that over space. Let us consider the case where the number of individuals does not change and without loss of generality consider the case with a single individual. Here we are considering situations as described in Section II, for example, within micro-fluidic devices where after each cell division, one daughter cell is removed from the ensemble under observation. Now we may consider the probability of making an observation within the interval : . Let us make observations on a single individual, and calculate the measured probability density for the interval : , where is the number of observations which fell within . These calculated probability densities follow the binomial distribution: . The root square difference between and the mean, , is .
We may readily construct the distribution defined within each interval as above. For suitably well behaved for which as , we may conclude that for any there exists a such that the root square difference between and tends to some limit below as . Smoothness of is more than sufficient and general enough for our discussion, so we will drop the notation and conclude that for smooth :
Thus we see that the ensemble behaves the same way over time and space: if we observe a single individual at many times, it recapitulates the behavior of the entire ensemble at a single time. We have shown that Eq. 14 is implied by Eq. 13; however, we have only shown this is true for the case where the number of individuals in the ensemble does not vary. We will show that in general, neither Eq. 14 nor Eq. 13 is implied by the other and will refer to them as follows:
Ergodicity Type I
We will define an ensemble to be of ergodicity type I if:
In other words, the ensemble has an unbiased selection of states which does not depend on time or the individual observed. Our condition of interest, is equivalent to Eq. 15
Ergodicity Type II
We will define an ensemble to be of ergodicity type II if:
as where , defined above, is the estimation of made from observations of any single individual. In other words, each individual in the ensemble will, over time, recapitulate the behavior of the entire ensemble.
We have just seen when the number of individuals is fixed, ergodicity type I implies ergodicity type II; however the converse is not true. Consider an ensemble composed of two individuals and which may occupy two possible states and . Suppose that at odd observations occupies and occupies and vice versa. The probability of observing either state in the entire ensemble is at any time and for both individuals; however for all times. Thus this ensemble is ergodic in the second sense but not ergodic in the first sense. See Fig. 3 for a cartoon illustrating the two types of ergodicity discussed.
When the number of individuals varies, as it does in the exponentially growing ensemble, ergodicity I does not imply ergodicity II. Let us turn to our ensemble of interest, the cell cycle duration distribution. We may utilize the von Foerster equation as we did earlier and note that using the same formalism assumes ergodicity I through the condition as we discuss above. We may return to the expression for the growth rate, : and note that is a convex function. By Jensen’s inequality this implies . Further, is strictly convex, which means equality holds only when . Considering the case where :
Considering the case where is not a Dirac delta function but a distribution with the same mean:
Thus, for any non-delta function distribution, . This means that does not tend to as the mean of tends to a value below the mean of . Thus ergodicity type I does not imply ergodicity type II when the number of individuals varies. In fact, for this example, ergodicity type I and ergodicity type II are mutually exclusive. This brings us to a statement of the relative strengths of these conditions (Please note that for the case of the exponentially growing ensemble discussed above, many observations of a single individual refers to a lineage of cell cycle durations, for example, current duration, mother cell duration, grandmother cell duration,…):
Hierarchy of Ergodicities Types I and II
For an ensemble with a fixed number of individuals, ergodicity type I implies ergodicity type II; however, ergodicity type II does not not imply ergodicity type I. Thus for these ensembles, ergodicity type I is the stronger condition. For an ensemble with a variable number of individuals, being of ergodicity type I does not necessarily imply the ensemble is of ergodicity type II and vice versa. Furthermore, for some ensembles, including the cell cycle duration distribution discussed above, they are mutually exclusive conditions. This distinction is of interest relative to the qualitative notion that an ergodic ensemble, recapitulates within a single individual over many observations the behavior of the entire ensemble at a single observation, or behaves the same way over time and space. These statements represent ergodicity type II and may not accurately describe an ensemble of ergodicity type I with a varying number of individuals.
Iv Growth Rate Gain
We have discussed how the cell cycle duration distribution is of ergodicity I under the condition which is assumed in the traditional formalism used to solve the von Foerster equation. Thus we know unless , the growth rate is higher than the “naive” growth rate and the bulk doubling time is shorter than the mean division time. This phenomenon is commonly referred to as the “growth rate gain”. It was discussed first as far back as the 1950’spowell1956growth and has been the object of renewed recent interestcerulus2016noise; hashimoto2016noise. At the center of the issue sits the finding that when the duration of mother-daughter cell cycles are positively correlated, the growth rate increases and when they are negatively correlated, the growth rate decreases. We have shown that even in the case without explicit correlation, ergodicity type I, the growth rate gain appears. This is due to the nature of the growing ensemble. When cells which divide quickly are equally likely to form daughter cells which divide quickly as they are to form daughter cells which divide slowly, this leads to an unequally large number of short cell cycles represented in the ensemble (see Fig. 3). On the other hand, we will show if the the ensemble is of ergodicity II, then there is no growth rate gain: the bulk doubling time is equal to the mean cell cycle duration.
Consider an exponentially growing ensemble which begins with a single individual such that all individuals in the ensemble will be of ergodicity II. Suppose at any given time, the ensemble contains cell cycle durations structured into lineages each of length . Let us consider the “error” of a single lineage, the root square difference between the mean of the lineage and , to be . The composite error of all the lineages within the ensemble labeling , and is:
Thus if all lineages are of ergodicity type II, the ensemble must also be of ergodicity type II. The bulk doubling time is is the mean of every cell cycle observed in the population (over some window of time observed) - which if the ensemble is of ergodicity II, is simply .
V Explicit Correlation
We have shown when ergodicity I is maintained in the case of the cell cycle duration distribution, , the ensemble growth rate is higher than . This growth rate gain may be attributed to the prevalence of lineages largely comprising cell cycles shorter than the mean. While there is no explicit mother-daughter correlation within the selection of cell cycle durations, there is an effective correlation arising from this variable weighting of lineages. We sought to probe the degree of this effective correlation through the addition of explicitly negative mother-daughter correlation. In this way, we could find the degree of negative correlation which must be imposed to return an ensemble to the naive growth rate, . We chose a form for based on the model presented inrochman2016grow:
The autocorrelation function generated from is:
The object of interest here is , the mother-daughter correlation. We simulated lineage trees of cell cycle durations with each successive cycle duration drawn from . See Fig. 4. The first generation cycle duration was drawn from a Gaussian distribution with mean and the standard deviation corresponding to the stationary distribution. Without loss of generality, was taken to be (unitless) and and were chosen to obtain the desired correlation and . Each cell lineage simulation was continued until a generation was reached where the cell cycle duration distribution was sufficiently close to the steady-state distribution evolving from . The test distribution was considered to be sufficiently close to the steady-state when the relative root square difference was no more than . Once the steady-state was reached, the growth rate of the population was calculated as the slope of the best-fit line to the logarithm of average cell number over trials.
We found a very substantial degree of negative correlation ( must be imposed to return an ensemble to the naive growth rate, ; however, this phenomenon depends on the of the cell cycle duration distribution since the growth rate gain is higher when the is larger. In other words, for an ensemble to attain a bulk growth rate of given only mother-daughter correlations (i.e. grandmother-granddaughter cell cycle durations are uncorrelated), an anti-correlation of about must be imposed. Under these conditions, a cell which divided after a cycle of duration is most likely to form daughter cells which attain cycles of duration . Stated another way, this implies that ensembles which enforce ergodicity I are essentially approximately correlated.
We have seen that the concept of ergodicity bifurcates into two distinct properties within ensembles that have a variable number of individuals: ergodicity I, the unbiased selection of states and ergodicity II, the recapitulation of whole-ensemble behavior from any single individual. Many of the nice properties of the traditional formalism used to solve the von Foerster equation rely on the assertion that the cell cycle duration distribution is of ergodicity I. Under these conditions, the bulk growth rate of the population is higher than due to a disproportionately large number of short cell cycles represented in the ensemble. When the ensemble is of ergodicity II, similar to the condition where there is a negative mother-daughter cell cycle duration correlation, the growth rate returns to . This observation corroborates the idea that the relationship between the variability of the cell cycle duration distribution and the population growth rate is impacted by mother-daughter correlationslin2017effects. Since the “growth rate gain” observed within ensembles of ergodicity I stems from disparities between lineages composed of primarily short cell cycle durations and those of long cell cycle durations, the larger the of the cell cycle duration distribution, , the greater the effect. In the most basic sense, this growth rate gain comes from asymmetric divisions. When ergodicity II is not enforced, lineages of exclusively short cell cycles arise. Many sister cells of cells within these lineages have long cell cycle durations. This phenomenon may be purely due to stochasticity present after division has finished or it might be due to programmatic asymmetry in the allocation of resources to each daughter cell. The use of the “FUCCI” live-cell cell cycle phase labeling systemsakaue2008visualizing will continue to provide the opportunity for experimental validation of these concepts relating to the age structure of the populationsandler2015lineage. In particular it may be interesting to evaluate how the age distribution of a monolayer varies as it closes a wound. In this case, one would expect cells on the edge of the wound to assume a different age distribution as they divide more frequently or migrate faster than cells far from the wound. Recently, related perspectives on how ergodicity and the lineage structure impacts growth rate have brought more interest to this topicthomas2017making; thomas2017single.
Asymmetric division, often a complex processmena2017asymmetric, has been well established in a variety of pro- and eukaryotesli2013art including the orchestration of stem cell differentiation and self-renewal. It has been shown that even E. coli display complex polar protein localizationlybarger2001polarity and that pathological polar aggregates can be asymmetrically inherited which may increase fitness by “rejuvenating” the daughter cell that accepts less damagewinkler2010quantitative. Completely symmetric division requires cells to fix inherited damage; otherwise, all cells will eventually have accumulated critical amounts. There are likely to be costs associated with coordinating asymmetric division and inherited damage mitigation which are balanced in optimal growth strategies. It has been reported that under some conditions, E. coli age within ten generationsstewart2005aging while under others, stable growth has been observed for hundreds of generationswang2010robust. Perhaps in the latter case, higher growth rates can be achieved through damage mitigation in all cells than the asymmetric inheritance of damage and subsequent loss of the damaged population.
NDR conducted the analysis in sections II, III, and IV. DMP completed the simulations in section V. NDR, DMP, and SXS wrote the paper.