Statistical Theory of Asymmetric Damage Segregation in Clonal Cell Populations

Many unicellular organisms among bacteria and yeasts proliferate by binary fission so that each mother cells divides into two seemingly identical daughter cells. However, more careful examination of cell lineages [1, 2] revealed that the progeny is in fact slightly asymmetric, one daughter cell has a slightly longer lifespan than the other. This difference, at least in rode-shaped bacteria E. coli, is apparently rooted in the fact that every bacterium is itself slightly asymmetric since its two poles always have a different age: the old pole that existed in its mother cell and the new pole that was created when its mother cell divided. The “age” of a cell can be defined as the age of the old pole (the number of generations it has been in existence). It was found that the daughter cell that inherits the older pole from its mother, grows slower and divides slightly later than the daughter cell that inherits the newer pole. A plausible hypothesis that can explain such a correlation is that the cells gradually accumulate damaged proteins which aggregates near their poles, and therefore the daughter cell that inherits the older pole, will also inherit a larger fraction of damaged proteins hitherto accumulated in the mother cell. Direct visualization of protein aggregates in growing cell lineages corroborates this conjecture [3, 4, 5, 6].

A number of computational and theoretical studies addressed the dynamics of asymmetric damage segregation (ADS) in a growing microbial population and their implications for the overall population fitness [7, 8, 9, 10, 11, 12, 13, 14, 15]. In [9, 11] it was demonstrated numerically that asymmetry accelerates the average growth of the population as a whole. In [14], a kinetic description of damage segregation was developed on the basis of a transport equation for the time-dependent damage distribution function that was applied both to deterministic and stochastic damage synthesis and segregation. Using that equation, moments and the average population growth rate were computed analytically in the limit of small asymmetry. However, a comprehensive understanding of the structure of damage distribution in a population of asymmetrically dividing cells is still lacking. In this paper, we focus on deterministic asymmetric division and analyze this structure using a general Frobenius-Perron equation for the at-birth damage distribution function. It turns out that if the rules controlling damage accumulation and inheritance are deterministic, the system governing the damage distribution is analogous, and for certain class of linear damage accumulation and division rules exactly equivalent, to the Iterated Function System (IFS) known in the theory of fractals [16]. Exploiting this analogy, we show that the asymptotic stationary distribution of damage is indeed fractal and find the spectrum of its fractal dimensions. For more general nonlinear models of damage accumulations and segregation we analyze the structure of the damage distribution functions numerically and show that the its fractal nature is robust. We discuss also deterministic irregularity of the damages in a lineage and relate it, for the linear damage redistribution rule, to chaotic properties of the IFS. If the damage accumulation and/or segregation have stochastic components, the distribution smears out but remains multi-peaked.

In this section we introduce a class of models for the asymmetric damage segregation that will form the basis of our theory. We will suppose that a cell is created at time $t_{0}$ and divides at time $t_{0}+\tau$. We denote the instantaneous cell damage (it can be, for example, the amount of damaged proteins, which is assumed to be a real number) $D(t)$, where $t_{0}\leq t\leq t_{0}+\tau$. The initial cell damage (inherited from the mother cell) is $x=D(t_{0})$, and the final damage just before division is $y=D(t_{0}+\tau)$, where $\tau$ is the cell lifetime (the latter is constant in some setups or damage-dependent in other formulations).

There are three components in each model that determines the distribution of damage within the population of growing and dividing cells:

1. 1.

   Damage gain specifies temporal dynamics of damage $D(t;x)$ in a cell starting from its value $x$ at birth at $t=t_{0}$ to division at $t=t_{0}+\tau$. We will denote the final cell damage as $y=D(t_{0}+\tau;x)$.

   Examples:

   1. a.

      Vedel et. al model [11]. In this model it it assumed that every cell adds a fixed amount of damage $\lambda$ over its lifetime, $y=x+\lambda$. The damage is added with the constant rate $\gamma$, so that for a lifetime $\tau$ this rate is $\gamma=\lambda/\tau$. The evolution of the damage is thus $D(t;x)=x+\lambda(t-t_{0})/\tau$.

   2. b.

      Damage obeys a linear differential equation $\dot{D}=\beta D+\gamma$, where the term $\beta D$ describes auto-catalytic production or degradation of damage (depending on the sign of $\beta$), with the initial condition $D(t_{0})=x$.

2. 2.

   Lifetime $\tau$ might be constant (damage-independent), or it can be a functional of the damage $D(t)$. Since in deterministic models at any moment of time the current damage $D(t;x)$ is pre-determined by the initial value $x$ at the instant of birth $t_{0}$, we can assume that the lifetime is a function of the initial damage, $\tau(x)$.

   Examples:

   1. a.

      In model [11], the lifetime is assumed to be a linear function of the cell damage at the division $y$. Because $y$ is determined by the damage at birth $x$, one can generally write $\tau(x)=\tau_{0}+\mu y(x)$.

   2. b.

      Chao model [9]. Here it is assumed that the lifetime is the time when some product $P$, whose synthesis rate depends on the damage, $\dot{P}=1-sD(t;x)$ with $P(t_{0})=0$, reaches a threshold value $P_{0}$. This leads to a nontrivial dependence of the cell lifetime on its initial damage $\tau(x)$, see Eq. (37) below.

3. 3.

   Damage inheritance: A deterministic rule according to which the damage of a mother cell is distributed between two daughters (no damage losses), $x_{1}=f(y),\ \ x_{2}=y-f(y),\quad 0\leq f(y)\leq y$.

   Examples:

   1. a.

      In model [9], a linear mapping $x_{1,2}=\frac{1\pm a}{2}y$ with a constant asymmetry parameter $0<a<1$ is adopted.

   2. b.

      In model [11], a nonlinear mapping $x_{1,2}=\frac{1\pm a(y)}{2}y$ with the asymmetry parameter being a function of the damage of the mother cell, $a(y)$, is suggested.

All model variants listed above are deterministic, however we will generalize our description to account for stochasticity in ADS in Sec. 8.

It is straightforward to simulate any specific ADS model numerically, starting, for example, from a single cell with a certain initial damage. Fig. 1a,b illustrate these dynamics for a model with the fixed damage gain (rule 1a of the previous Section), linear lifetime function (2a), and linear damage inheritance rule (3a). In the following we will label this model as 1a2a3a. Fig. 1a corresponds to the constant lifetime ($\mu=0$) and Fig. 1b to the damage-dependent lifetime ($\mu>0$). Note that since the number of cells in such a simulation grows exponentially with time, it cannot be extended very far. As an alternative, one can simulate a Moran-type process with a constant population of $N$ cells when at each cell division, another random cell is removed from the population, so the population size remains constant. The statistical features of both processes are the same in the “thermodynamic limit” $N\to\infty$, where one can neglect finite-size fluctuations. We can also expect them to be similar for finite $N$ at the time when the size of a growing population is also $N$.

The main difference between the cases in panels (a) and (b) is that for the damage-independent lifetime all cells divide simultaneously, while for damage-dependent lifetimes, synchrony of division times is quickly lost. In both cases, after a short transient, the distribution of damage in just-born cells appears to reach a broad and multi-peaked stationary state (Fig. 1c,d) with a finite support between certain $x_{min}$ and $x_{max}$. It is straightforward to find $x_{min}$ and $x_{max}$ for a given model as the asymptotic limit of inheriting only smaller or only larger fractions of mother’s cell damage within a lineage. For example, for Model 1a2a3a, $x_{min}=\lambda\frac{1-a}{1+a}$, $x_{max}=\lambda\frac{1+a}{1-a}$.

It is, however, not trivial to interpret these simulations. Indeed, the process is fully deterministic: an initial damage of a cell fully determines damages and fission times of all its descendants. But it is not a classical deterministic dynamical system because the number of cells grows, and we cannot represent the time evolution of the damage values in all descendant cells as a single trajectory. However, while each mother cell has two daughters, the daughter cell damage uniquely specifies the damage of its mother, at least for the model 1a2a3a, if parameters $a$ and $\lambda$ are also known. The reason for this is that the damages of daughters receiving larger and smaller fractions of the mother’s damage belong to non-overlapping sets $(x_{min},\lambda)$ and $(\lambda,x_{max})$ (this is not the case for all ADS models, see below). Thus, for the model 1a2a3a we can construct a deterministic ‘back-in-time’ map for the damage $x(n)$ of the just-born cells in $n$-th generation

(Fig. 2a) and generate a unique “lineage” damage trajectory back in time, as illustrated by Fig. 2b where we concatenated the damage time courses in different generations.

The main conclusion of Fig. 2b is that the deterministic nature of the model notwithstanding, the time course of the damage is highly irregular, which suggests that this process may be approached from the viewpoint of the chaos theory. The mapping (1) is indeed an expanding piecewise linear map which is a prototypical model of deterministic chaos [17]. This immediately explains the irregularity of the lineage damage trajectory. In fact, relation (1) can be considered a Poincaré map for the continuous trajectories of the full time-dependent damage dynamics. The latter is described not by differential equations, but rather by a combination of a continuous evolution of damage during the cell lifetime interrupted by a discrete transition from mother to daughter by the rule of damage inheritance. In this sense our hybrid continuous-discrete system resembles an integrate-and-fire model for firing neurons, where a continuous voltage increase (current integration) is combined with a jump rule (firing of a spike and a reset of the voltage).

We stress here that although a back-in-time trajectory can be always extracted from a forward-in-time simulation of a population of dividing cells, not in every case this back-in-time trajectory can be associated with a deterministic dynamical system like the map (1). It requires the uniqueness of the back-in-time map, i.e. for any initial damage of a daughter cell there should be only one possible value of the mother’s initial damage. A detailed analysis (Sec. 5) shows that this is not always the case. For example, for model 1b2a3a with $\beta>0$, two branches of the back-in-time map similar to Fig. 2b) overlap, and so neither forward-in-time, nor back-in-time lineage trajectories can be computed by iterating a deterministic 1D map.

The essential ingredient of any ADS model is the dependence of cell lifetime on its damage. Continuing with our analogy between the mother-to-daughter damage inheritance and a Poincaré map, one can interpret the lifetimes as Poincaré return times. The properties of these times affect the regularity of the damage trajectories Fig. 2a; the corresponding notion in the chaos theory is phase coherence [18]. One may distinguish a generic case of damage-dependent lifetimes from a degenerate case of fixed lifetimes. For fixed lifetimes, the time evolution of the damage is only partially irregular: Divisions occur at regular time intervals (Fig. 1a), and there is no mixing for the continuous-time dynamical process. All lineages starting from the same cell undergo divisions simultaneously, and so the mean population damage continues to oscillate indefinitely (see Fig. 3a, red line). In other words, all the cells remain “in phase”. In contrast, for a generic damage dependence of lifetimes, the intervals between sequential divisions are chaotic (since the initial damages of cells in a lineage are chaotic), and the continuous-time process is mixing (Fig. 1b). Thus, different lineages starting from the came cell decorrelate (in other words, the phases of different cells become scattered) and the mean population damage eventually settles into a steady state (see Fig. 3a, blue line). Another way to quantify this difference is to compute the normalized autocorrelation function of the damage trajectory

where the angular brackets denote averaging over different lineages ¹¹1 For averaging we used here so-called retrospective (or historic) sampling by taking lineages of all cells present at a given time (see [19, 20]). For an alternative chronological (or forward) sampling), when each daughter of a given cell is selected with equal probability, the averages may be slightly different since it overestimates older, slower-dividing cells.

In Fig. 3b we compare this autocorrelation function for model 1a2a3a with $\mu\neq 0$ and $\mu=0$. One can see that while for damage-dependent lifetimes ($\mu\neq 0$) this function decays to zero, for constant lifetimes it initially decays (reflecting the irregularity of the damages in just-born cells), but asymptotically at large time lags it oscillates periodically without decay, reflecting phase coherence of the corresponding damage trajectory. For the latter case, the autocorrelation function can be computed analytically (see A).

In the previous section we demonstrated that the dynamics of the asymmetric damage segregation in cell lineages is irregular, thus it appears appropriate to describe them statistically. In this section, we formulate a kinetic description of ADS in the thermodynamic limit of a large number of cells. It is similar to the recently developed transport approach [15, 14] but differs in some details.

As we have argued above, in many aspects the damage dynamics is analogous to the dynamics of chaotic oscillators, and we will follow this analogy in the construction of the statistical theory. It is convenient to introduce a phase of the cell cycle $\phi$ that changes linearly from $\phi=0$ (birth) to $\phi=1$ (division) with the rate $\omega$ (which is simply the inverse of the cell lifetime $\tau$) that generally may depend on the initial cell damage $x$:

Next, instead of $D(t;x)$ above, we can introduce the damage during the cell lifetime as a deterministic function of the initial state $x$ and the phase $\phi$, $D=F(x,\phi)$. This function is related to the function $D(t;x)$ above as

Let us now introduce a two-dimensional population density $N(x,\phi,t)$, so that $N(x,\phi,t)dxd\phi$ is the number of cells in a small interval of the initial damages $dx$ and in a small interval of phases $d\phi$. With a proper normalization, the quantity $N(x,\phi,t)$ can be interpreted as a probability density in the space of deterministic variables $x,\phi$. We stress here that the irregularity of the dynamics is not needed for this: even regular deterministic processes (e.g. populations of periodic oscillators) can be described via a probability density in the phase space, with an additional assumption of random initial phases in an ensemble of copies [21]. This quantity is similar to the distribution $\psi(u,x,t)$ introduced in [14], which instead of our phase $\phi$ depends on the actual cell age $u$ measured in units of time. We believe that our distribution of cell ages in terms of the phase $\phi$ (i.e. the fraction of the cell’s lifetime) has its advantages, especially when the lifetime varies strongly from cell to cell.

Instead of dealing with an exponentially growing population of dividing cells, it is convenient to introduce cell death, so that a certain fraction of cells $\Gamma N(x,\phi,t)dt$ dies within a small time interval $dt$ irrespective of $x$ or $\phi$. The death rate $\Gamma$ in principle can be arbitrary, but if we choose it exactly balance the (yet unknown) average asymptotic proliferation rate, the total number of cells will remain constant, and the distribution $N(x,\phi,t)$ should eventually approach a stationary state $N(x,\phi)$. It is easy to see that the death rate $\Gamma$ found from the stationarity condition of the asymptotic solution for a model with death is exactly the growth rate of a proliferating cell population without death. Alternatively, we can assume a Moran-type process [22] when at the time of every binary fission, one other randomly chosen cell is taken away, so the total number of cells remains precisely constant at all times, not just asymptotically and in the thermodynamic limit. It can be shown (see B) that asymptotically in the large-time limit the kinetic equation for the Moran-type process is equivalent to the one presented here.

The transport equation for $N(x,\phi,t)$ in the system with cell death follows from the probability conservation:

where the second term on the l.h.s. describes cells advection from phase 0 to phase 1 with phase velocity $\omega(x)$ according to (3), and the r.h.s. term describes cell death. This equation formally coincides with the transport equation for the normalized damage distribution in an exponentially growing population [14] with $\Gamma$ playing the role of the population growth rate instead of the death rate in our model. We assume that the death rate $\Gamma$ is indeed chosen properly to balance the population growth. In this case, and if the initial phases in a population are random, the probability distribution asymptotically approaches a stationary state described by the equation ²²2If the dynamics are mixing, the initial random phase approximation is not required, and the stationary asymptotic distribution will be achieved from an arbitrary initial distribution, but for the degenerate case of damage-independent lifetimes, initial “randomization” is necessary.

The solution of this equation describes the exponentially decaying flux of cells from $\phi=0$ towards $\phi=1$:

To close the system and to ensure stationarity, we need to add the boundary (transition) condition balancing the incoming flux of the daughter cells $N(\bar{x},0)\omega(\bar{x})$ at $\phi=0$ with the outgoing flux of mother cells $N(x,1)\omega(x)$ at $\phi=1$ (due to the cell division, the number of cells leaving the interval at $\phi=1$ has to be exactly the half of the number of cells entering at $\phi=0$, thus factor 2 in Eq. (14) below) .

The connection is determined by the rule specifying how the damage at the end of the life cycle $y=F(x,1)$ is redistributed within daughter cells. For a general damage inheritance rule (see 3. in section 2) we have

For monotonous functions $F$ and $f$, depending on their specific forms, a daughter cell with initial damage $\bar{x}$ can descend from one of two mothers with initial damages $x_{1,2}=F^{-1}(f_{1,2}^{-1}(\bar{x}),1)$.

The balance condition for the fluxes thus can be written as

Dividing both sides by $dx$ we arrive at the basic Frobenius-Perron operator for the damage distribution in a population of dividing cells at the beginning of their cell cycle

Alternatively, we can write this Frobenius-Perron equation for the product $P(x)=\omega(x)N(x,0)$ in the following integral form

where $g_{1,2}(x)=f_{1,2}(F(x,1))$ and we dropped the bar over $x$. It is easy to see that properly normalized $P(x)$ is the density of just-born cells (the fraction of the whole population of cells that that are born within a small unit time interval with initial damage $x$).

This equation can also be deduced directly from the general equation for the stationary initial damage distribution function derived in [14]

where $f(x|x^{\prime})$ is the probability to have initial damage $x$ in a cell if the initial damage in its mother cell is $x^{\prime}$. Eq. (12) can describe both deterministic and stochastic damage distribution scenarios. To obtain our deterministic (11) one just needs to take $g(x|x^{\prime})={1\over 2}[\delta(x-g_{1}(x^{\prime}))+\delta(x-g_{2}(x^{\prime}))]$.

The conservation of the number of cells requires

or

This self-consistency condition together with Eq. (11) uniquely determines the death/growth rate $\Gamma$.

Asymptotically, the solution of the Frobenius-Perron equation (10) (provided the death rate $\Gamma$ satisfies (13)) converges to the stationary damage distribution at the moments of birth $N(x,0)$. It is easy to see that this distribution has a finite support between $x_{min}$ and $x_{max}$ which are determined from equations $x_{min}=g_{1}(x_{min}),\ x_{max}=g_{2}(x_{max})$ (we assume for definiteness that $g_{1}(x)<g_{2}(x)$).

Let us introduce the total number of cells (here we use also (7)):

Then the stationary probability distribution density over initial damages and phases is

We can use this two-dimensional density to obtain the density for the distribution of the damages (now not initial one, but just observed at all possible phases), by the following expression

Here we used expression (4) which relates the damage at the phase $\phi$ to the initial damage $x$. We can also obtain he distribution over the phases $w(\phi)$ as the marginal distribution by integrating the two-dimensional density (16)

In this section we use the kinetic theory developed in the previous section to analyze the statistical properties of the model with a fixed lifetime $\tau=T$ and linear mappings

with constant $A$ and $B$ where for definiteness, we assume $A<B$.

It is easy to see that linear transformations (18) with $A+B=1$ correspond to rules 1a and 3a of Sec. 2, where

The same linear mappings also can be deduced from rule 1b, however in this case, integration of the ODE for the damage within the cell cycle yields

and so for non-zero $\beta$, $A+B\neq 1$. If $\beta<0$, then $A+B<1$, and if $\beta>0$, then $A+B>1$. Note that for the damage to remain bounded at all times it is required that $A<1,B<1$, i.e. $A+B<2$.

The fixed lifetime may correspond either to rule 2a with $\mu=0$ (then $T=\tau_{0}$) or to rule 2b with $s=0$ (then $T=P_{0}$).

Here, we consider the same linear model with the mother-daughter damage inheritance relations (18), but assume that cell lifetimes depend on the damage, $\tau(x)$.