A transport approach to relate asymmetric protein segregation and population growth

Genetically identical unicellular organisms may nevertheless exhibit different traits such as growth rate or generation time, a phenomenon known as phenotypic variability [1]. One cause of such diversity is the asymmetric allocation of key proteins at division. For instance, some bacteria and eukaryotic cells accumulate protein aggregates when exposed to stress and tend to segregate more of them to one of the two newborn cells at a division event. As a result, the 'rejuvenated' cells (i.e. those with less damage) grow faster than the 'senescent' cells [2–7]. There have been a number of theoretical works [8–11] trying to quantify the effect of asymmetric segregation of proteins to the overall growth of the population. In this paper, we strengthen the theory developed by Lin et al [8] to relate the asymmetric protein segregation and the population growth rate by implementing a more general method involving transport equations and moment expansions. This allows us to recover previous results rigorously, and go beyond them. We extend the model to consider stochastic asymmetric segregation. Lastly, we compare alternative models of asymmetric segregation by applying the self-consistent equations.

We define a generalization of the model proposed by Lin et al [8]. Consider a clonal microbial population that grows exponentially with rate Λ_p . We assume that as a cell grows by volume ΔV, the amount of key proteins D increases by ΔD = SΔV. Each cell's inter-division time (τ) is a function of the initial concentration of the key protein (σ_b). For simplicity, suppose every cell is born with volume V_b = 1 and divides when the volume reaches V_d = 2. As it divides, the key proteins segregate asymmetrically according to the asymmetry parameter a: one daughter cell inherits $\frac{1+a}{2}$ of the key proteins from the mother and the other gets the rest. We call the key protein damaging if the generation time τ[σ_b] is an increasing function of the initial protein concentration and beneficial if it is a decreasing function (figure 1).

Zoom In Zoom Out Reset image size

While our analysis can be applied to any analytical function τ[σ_b], we assume the instantaneous growth rate λ[σ] is a Hill-type function of the key protein concentration of the cell, as proposed by Lin et al [8]. These assumptions by Lin et al are based on observations made in experimental studies, where the growth rate function of microbes often has an inflection point [5, 6].

$$\begin{equation}\lambda [\sigma ]=\frac{{\lambda }_{0}+{\lambda }_{1}{\sigma }^{n}}{1+{\sigma }^{n}},\end{equation} \tag{ 1 }$$

$$\begin{align}\tau [{\sigma }_{\mathrm{b}}]& ={\int }_{1}^{2}\frac{1}{V\lambda [\sigma ]}\enspace \mathrm{d}V\\ & ={\int }_{1}^{2}\frac{1}{V\lambda [({\sigma }_{\mathrm{b}}-S)/V+S]}\enspace \mathrm{d}V.\end{align} \tag{ 2 }$$

Every time a cell divides, the asymmetry parameter a is randomly and independently drawn from a probability distribution g(a). The deterministic model by Lin et al [8] assumes a delta function for g(a), but here we allow a to fluctuate.

3.1. Transport equation

Recently, Levien et al [12] showed that a transport equation can be used to get a self-consistent equation for the population distribution of phenotype and population growth rate even when the properties of the cells are correlated across the generations. Here we follow their derivation by setting the phenotype of interest as the key protein concentration at birth σ_b.

Assume that the population size is N at time t and it grows exponentially with the rate Λ_p . Let us also assume that the population has reached a stable probability distribution ψ(u, σ_b, t), where u is time since birth, σ_b is the key protein concentration at birth, and t is the global time. The probability distribution then evolves during an infinitesimal time step dt as:

$$\begin{align}N\psi (u,{\sigma }_{\mathrm{b}},t)& =N\enspace {\mathrm{e}}^{-{{\Lambda}}_{p}\mathrm{d}t}\psi (u-\mathrm{d}t,{\sigma }_{\mathrm{b}},t-\mathrm{d}t),\\ & =N(1-{{\Lambda}}_{p}\enspace \mathrm{d}t)\left(\psi (u,{\sigma }_{\mathrm{b}},t)-\mathrm{d}t\enspace \frac{\partial \psi }{\partial u}-\mathrm{d}t\enspace \frac{\partial \psi }{\partial t}\right).\end{align} \tag{ 3 }$$

Rearranging the terms and dividing both sides by dt, we get

$$\begin{equation}\frac{\partial \psi }{\partial u}+\frac{\partial \psi }{\partial t}=-{{\Lambda}}_{p}\psi .\end{equation} \tag{ 4 }$$

In a steady state population, there is no dependence on t, so ∂_t ψ = 0. Therefore,

$$\begin{equation}\frac{\partial \psi (u,{\sigma }_{\mathrm{b}})}{\partial u}=-{{\Lambda}}_{p}\psi (u,{\sigma }_{\mathrm{b}}),\end{equation} \tag{ 5 }$$

$$\begin{equation}\psi (u,{\sigma }_{\mathrm{b}})=\psi (0,{\sigma }_{\mathrm{b}}){\mathrm{e}}^{-{{\Lambda}}_{p}u}.\end{equation} \tag{ 6 }$$

The boundary conditions are set by the division. The probability distribution of having a newborn cell with a protein concentration σ_b is

$$\begin{equation}\psi (0,{\sigma }_{\mathrm{b}})=2\iint f({\sigma }_{\mathrm{b}}\vert {\sigma }_{\mathrm{b}}^{\prime },{a}^{\prime })g({a}^{\prime })\psi (\tau [{\sigma }_{\mathrm{b}}^{\prime }],{\sigma }_{\mathrm{b}}^{\prime })\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }\enspace \mathrm{d}{a}^{\prime }.\end{equation} \tag{ 7 }$$

The probability distribution f(σ_b|σ_b', a') and g(a') depend on the physiological model we choose. For instance, in a deterministic model, we have

$$\begin{equation}f({\sigma }_{\mathrm{b}}\vert {\sigma }_{\mathrm{b}}^{\prime },{a}^{\prime })=\frac{1}{2}\delta \left({\sigma }_{\mathrm{b}}-\frac{1+{a}^{\prime }}{2}({\sigma }_{\mathrm{b}}^{\prime }+S)\right)+\frac{1}{2}\delta \left({\sigma }_{\mathrm{b}}^{\prime }-\frac{1-{a}^{\prime }}{2}({\sigma }_{\mathrm{b}}^{\prime }+S)\right),\end{equation} \tag{ 8 }$$

$$\begin{equation}g({a}^{\prime })=\delta ({a}^{\prime }-\overline{a}).\end{equation} \tag{ 9 }$$

Applying the solution given by equation (6) (in the general case), we can eliminate the dependence on u. Thus, we can rewrite to obtain a self-consistent equation:

$$\begin{equation}{\psi }_{\text{birth}}[{\sigma }_{\mathrm{b}}]=2\iint f({\sigma }_{\mathrm{b}}\vert {\sigma }_{\mathrm{b}}^{\prime },{a}^{\prime }){\psi }_{\text{birth}}[{\sigma }_{\mathrm{b}}^{\prime }]g({a}^{\prime }){\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\sigma }_{\mathrm{b}}^{\prime }]}\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }\enspace \mathrm{d}{a}^{\prime }.\end{equation} \tag{ 10 }$$

In the next section, we apply the self-consistent equation (10) to a toy model of E. coli recently published by Blitvic and Fernandez [13] as a simple example.

3.2. Toy model of the aging of E. coli [13]

Let us introduce a simple model published by Blitvic et al to demonstrate how to use our approach to connect the asymmetry between two newborn cells and population growth rate [13]. Blitvic et al assume that as an E. coli cell divides, one daughter cell is rejuvenated and the other senescent where the generation time of a rejuvenated cell (τ[R]) is shorter than a senescent cell (τ[S]). In this model, there are two phenotype x = S for a senescent cell and x = R for a rejuvenated one.

Then the transition density f(x|x') from a phenotype x' to x is

$$\begin{equation}f(x\vert {x}^{\prime })=\frac{1}{2}\delta (x-R)+\frac{1}{2}\delta (x-S).\end{equation} \tag{ 11 }$$

The self-consistent equation is

$$\begin{equation}\psi [x]=\int \psi [{x}^{\prime }]f(x\vert {x}^{\prime }){\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{x}^{\prime }]}\mathrm{d}{x}^{\prime }.\end{equation} \tag{ 12 }$$

Integrating both sides with respect to x, we have

$$\begin{equation}1={\mathrm{e}}^{-{{\Lambda}}_{p}\tau [R]}+{\mathrm{e}}^{-{{\Lambda}}_{p}\tau [S]}.\end{equation} \tag{ 13 }$$

Because the exponential function is convex, by Jensen's inequality,

$$\begin{equation}1={\mathrm{e}}^{-{{\Lambda}}_{p}\tau [S]}+{\mathrm{e}}^{-{{\Lambda}}_{p}\tau [R]}{ >}2\enspace {\mathrm{e}}^{-{{\Lambda}}_{p}\left\langle \tau [x]\right\rangle },\end{equation} \tag{ 14 }$$

where $\left\langle \tau [x]\right\rangle =\frac{1}{2}(\tau [S]+\tau [R])$. Thus, the asymmetry of generation time of two distinct phenotypes enhances the population growth rate if the mean of the generation time, $\left\langle \tau [x]\right\rangle$ is kept constant. Note that equation (13) was also derived in [14] in the context of the asymmetric division of the budding yeast Saccharomyces cerevisiae.

In the next section, we solve equation (10) numerically and analytically more generally.

3.3. Numerically solving the self-consistent equation

The self-consistent equation (10) can be numerically solved by iterating the following procedure. First, we start with initial guesses for the solutions ${\psi }_{\text{birth}}^{(0)}[{\sigma }_{\mathrm{b}}]$ and ${{\Lambda}}_{p}^{(0)}$. We then update our guess for the distribution ψ_birth[σ_b] using equation (10). Next, we update Λ_p by using the Euler–Lotka equation [12], obtained by integrating equation (10) with respect to σ_b. This procedure is summarized in this set of equations:

$$\begin{equation}{~{\psi }}_{\text{birth}}^{(n)}[{\sigma }_{\mathrm{b}}]=2\int f({\sigma }_{\mathrm{b}}\vert {\sigma }_{\mathrm{b}}^{\prime },{a}^{\prime }){\psi }_{\text{birth}}^{(n-1)}[{\sigma }_{\mathrm{b}}^{\prime }]g({a}^{\prime }){\mathrm{e}}^{-{{\Lambda}}_{p}^{(n-1)}\tau [{\sigma }_{\mathrm{b}}^{\prime }]}\mathrm{d}{a}^{\prime }\enspace \mathrm{d}{\sigma }_{\mathrm{b}}^{\prime },\end{equation} \tag{ 15 }$$

$$\begin{equation}{\psi }_{\text{birth}}^{(n)}[{\sigma }_{\mathrm{b}}]=\frac{{~{\psi }}_{\text{birth}}^{(n)}[{\sigma }_{\mathrm{b}}]}{\int \mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }{~{\psi }}_{\text{birth}}^{(n)}[{\sigma }_{\mathrm{b}}^{\prime }]},\end{equation} \tag{ 16 }$$

$$\begin{equation}1=2\iint {\psi }_{\text{birth}}^{(n)}[{\sigma }_{\mathrm{b}}^{\prime }]g({a}^{\prime }){\mathrm{e}}^{-{{\Lambda}}_{p}^{(n)}\tau [{\sigma }_{\mathrm{b}}^{\prime }]}\mathrm{d}{a}^{\prime }\enspace \mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }.\end{equation} \tag{ 17 }$$

${\psi }_{\text{birth}}^{(n)}[{\sigma }_{\mathrm{b}}]$ and ${{\Lambda}}_{p}^{(n)}$ converge to constant values as n → ∞. Note that the solution is not sensitive to the choice of ${\psi }_{\text{birth}}^{(0)}$ and ${{\Lambda}}_{p}^{(0)}$. In appendix D, we discuss the shape of the distribution ψ_birth[σ_b].

In figure 2, we compare the numerical solution for equation (10) in the case of deterministic a (i.e. g(a') is a delta function) and the fitted exponential growth rate obtained by simulating the population directly (see appendix A for details of the simulation).

Zoom In Zoom Out Reset image size

Figure 2 confirms that equation (10) has sufficient information for calculating the population growth rate. The iterative procedure is more accurate, faster, and less memory-intensive than a brute force simulation which involves exponentially growing population size. However, the main shortcoming of the numerical solution is that it does not give any insight about how Λ_p depends on S and a. It is possible infer the relation by finding Λ_p for various values of S and a, but for every parameter value, we have to go through equations (15)–(17), which is highly inefficient. For those reasons, the moment expansion method we discuss in the following section is a good alternative way of calculating Λ_p for small a.

3.4. Moment expansion

Equation (10) can be approximately solved by expanding the moments of ψ_birth[σ_b]. For simplicity, we start with the deterministic case. As we will explain in the following sections, it is straightforward to generalize this moment expansion method to the stochastic case.

When the key protein segregation is perfectly symmetric (i.e. a = 0), every cell has a concentration of key protein of σ_b = S at cell birth. Therefore, for a small a, the distribution ψ_birth[σ_b] should be sharply peaked around S. Based on this, we are motivated to expand the moment of the distribution around σ_b = S by defining Δσ_b = σ_b − S. In order to calculate the nth moment, we multiply both sides of equation (10) by ${\Delta}{\sigma }_{\mathrm{b}}^{n}$ and integrate them with respect to σ_b. We use equations (8) and (9) to specify f(σ_b|σ_b', a') and g(a'). Invoking Newton's binomial theorem, we obtain:

$$\begin{equation}\langle {\Delta}{\sigma }_{\mathrm{b}}^{n}\rangle =\frac{1}{{2}^{n-1}}\sum\limits _{m=0}^{\lfloor n/2\rfloor }\left(\genfrac{}{}{0pt}{}{n}{2m}\right){a}^{2m}\left\langle {\Delta}{\sigma }_{\mathrm{b}}^{n-2m}{({\Delta}{\sigma }_{\mathrm{b}}+2S)}^{2m}\enspace {\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\Delta}{\sigma }_{\mathrm{b}}+S]}\right\rangle .\end{equation} \tag{ 18 }$$

Note that by convention, $\left(\genfrac{}{}{0pt}{}{n}{0}\right)=1$ and ⌊n/2⌋ is the greatest integer that is less than or equal to n/2. From the right-hand side of equation (18), we can infer that $\langle {\Delta}{\sigma }_{\mathrm{b}}^{n}\rangle =O({a}^{2\lfloor n/2\rfloor })$. This means that by solving a system of 2k + 1 equations (n = 0, ..., 2k), we will find Λ_p to O(a^2k ) accuracy. For each equation, we need to expand the exponential function ${\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\Delta}{\sigma }_{\mathrm{b}}+S]}$ around Δσ_b = 0 to order a^2⌊n/2⌋. To illustrate the process more clearly, let us calculate Λ_p to O(a²) accuracy.

If a = 0, every cell has σ_b = S when the population is stabilized, so ${{\Lambda}}_{p}=\frac{\mathrm{ln}\enspace 2}{\tau [S]}$. Due to symmetry, Λ_p [a] = Λ_p [−a]. Thus, to solve Λ_p to O(a²) is equivalent to finding a coefficient C₂ defined as below,

$$\begin{equation}{{\Lambda}}_{p}=\frac{\mathrm{ln}\enspace 2}{\tau [S]}(1+{C}_{2}{a}^{2}+O({a}^{4})).\end{equation} \tag{ 19 }$$

Taylor expanding τ[Δσ_b + S] around Δσ_b = 0 to O(a²) we find

$$\begin{equation}\tau [{\Delta}{\sigma }_{\mathrm{b}}+S]=\tau [S]+{\tau }^{\prime }[S]{\Delta}{\sigma }_{\mathrm{b}}+\frac{{\tau }^{{\prime\prime}}[S]}{2}{\Delta}{\sigma }_{\mathrm{b}}^{2}.\end{equation} \tag{ 20 }$$

Expanding ${\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\Delta}{\sigma }_{\mathrm{b}}+S]}$ to O(a²) leads to:

$$\begin{equation}{\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\Delta}{\sigma }_{\mathrm{b}}+S]}=\frac{1}{2}\left[1-\mathrm{ln}\enspace 2\left({C}_{2}{a}^{2}+\frac{{\tau }^{\prime }[S]}{\tau [S]}{\Delta}{\sigma }_{\mathrm{b}}+\left(\frac{{\tau }^{{\prime\prime}}[S]}{2\tau [S]}-\frac{\mathrm{ln}\enspace 2}{2}\frac{{\tau }^{\prime }{[S]}^{2}}{\tau {[S]}^{2}}\right){\Delta}{\sigma }_{\mathrm{b}}^{2}\right)\right].\end{equation} \tag{ 21 }$$

Applying these to equation (18) while setting n = 0, 1, or 2, we can set up a system of three linear equations with the unknowns being ${C}_{2},\left\langle {\Delta}{\sigma }_{\mathrm{b}}\right\rangle$ and $\left\langle {\Delta}{\sigma }_{\mathrm{b}}^{2}\right\rangle$. (Note that here τ[S] is a known function specified in the model.)

$$\begin{equation}1=1-\mathrm{ln}\enspace 2\left[{C}_{2}{a}^{2}+\frac{{\tau }^{\prime }[S]}{\tau [S]}\langle {\Delta}{\sigma }_{\mathrm{b}}\rangle +\frac{{\tau }^{{\prime\prime}}[S]}{2\tau [S]}\langle {\Delta}{\sigma }_{\mathrm{b}}^{2}\rangle -\frac{\mathrm{ln}\enspace 2}{2}\frac{{\tau }^{\prime }{[S]}^{2}}{\tau {[S]}^{2}}\langle {\Delta}{\sigma }_{\mathrm{b}}^{2}\rangle \right],\end{equation} \tag{ 22 }$$

$$\begin{equation}\langle {\Delta}{\sigma }_{\mathrm{b}}\rangle =\frac{\langle {\Delta}{\sigma }_{\mathrm{b}}\rangle }{2}-\frac{\mathrm{ln}\enspace 2}{2}\frac{{\tau }^{\prime }[S]}{\tau [S]}\langle {\Delta}{\sigma }_{\mathrm{b}}^{2}\rangle ,\end{equation} \tag{ 23 }$$

$$\begin{equation}\langle {\Delta}{\sigma }_{\mathrm{b}}^{2}\rangle =\frac{\langle {\Delta}{\sigma }_{\mathrm{b}}^{2}\rangle }{4}+{S}^{2}{a}^{2}.\end{equation} \tag{ 24 }$$

Solving these, we get

$$\begin{equation}\langle {\Delta}{\sigma }_{\mathrm{b}}\rangle =-\frac{4\enspace \mathrm{ln}\enspace 2}{3}\frac{{\tau }^{\prime }[S]}{\tau [S]}{S}^{2}{a}^{2},\end{equation} \tag{ 25 }$$

$$\begin{equation}\langle {\Delta}{\sigma }_{\mathrm{b}}^{2}\rangle =\frac{4}{3}{S}^{2}{a}^{2},\end{equation} \tag{ 26 }$$

$$\begin{equation}{C}_{2}=-\frac{2}{3}\frac{{\tau }^{{\prime\prime}}[S]}{\tau [S]}{S}^{2}+2\enspace \mathrm{ln}\enspace 2\frac{{\tau }^{\prime }{[S]}^{2}}{{\tau }^{2}[S]}{S}^{2}.\end{equation} \tag{ 27 }$$

Therefore, the population growth rate is

$$\begin{equation}{{\Lambda}}_{p}=\frac{\mathrm{ln}\enspace 2}{\tau [S]}\left(1+\left(-\frac{2}{3}\frac{{\tau }^{{\prime\prime}}[S]}{\tau [S]}+2\enspace \mathrm{ln}\enspace 2\frac{{\tau }^{\prime }{[S]}^{2}}{{\tau }^{2}[S]}\right){S}^{2}{a}^{2}\right).\end{equation} \tag{ 28 }$$

This is general to any choice of τ[σ_b], but in particular, if the generation time is described by equation (2), we can rewrite the population growth rate in terms of λ[S].

$$\begin{equation}{{\Lambda}}_{p}=\lambda [S]\left(1+\frac{{\lambda }^{{\prime\prime}}[S]}{4\enspace \mathrm{ln}\enspace 2\lambda [S]}{S}^{2}{a}^{2}\right).\end{equation} \tag{ 29 }$$

This second order solution for Λ_p is consistent with the result by Lin et al [8] but is derived here more rigorously and more generally since it applies to any functional form of τ[σ_b]. Furthermore, we can calculate Λ_p to an arbitrarily high order as desired by solving for ${{\Lambda}}_{p}=\frac{\mathrm{ln}\enspace 2}{\tau [S]}\left(1+{C}_{2}{a}^{2}+{C}_{4}{a}^{4}+\cdots \enspace \right)$, as we elaborate in appendix C.

In the previous sections, we explained how to solve equation (10) and predict the population growth rate Λ_p as a function of growth rate of a single cell λ[σ], protein accumulation rate S, and asymmetry a. In the context of evolution, we can imagine that a single-celled organism may adapt to environmental stress (characterized by S and λ[σ] in our model) by tuning the asymmetry a. In population genetics, fitness of a population is often defined as an exponential growth rate of the population size, and traits that increase the fitness tend to survive longer due to natural selection [15]. In fact, there are several experimental studies reporting that segregation of protein aggregates tend to become more asymmetric as temperature increases above the optimal temperature [2–5]. In this section, we present how to find the optimal asymmetry a_c, the value of a that maximizes the population growth rate Λ_p , as a function of population accumulation rate S. As Lin et al [8] did, we find a sharp phase transition of a_c if the key protein is deleterious and a smooth transition in the case of beneficial key protein.

Previously, Lin et al [8] inferred a_c[S] by interpolating the solutions near a = 0 and at a = 1 with a minimal number of extrema. Lin et al used a Landau approach to find the population growth rate Λ_p to O(a²) accuracy for a small a. In summary, the Landau approach investigates the shape of the fitness curve f = Λ_p [a] − Λ_p [a = 0] as a function of a at a given key protein accumulation rate S (i.e. the fitness is analogous to a free energy with an order parameter a, and S to temperature). They derived the relationship between the fitness f and asymmetry a from an equation for the total cell volume, V = ∑_i v_i (t).

$$\begin{equation}\frac{\mathrm{d}V}{\mathrm{d}t}=\sum\limits _{i}\lambda [{\sigma }_{\mathrm{b},i}]{v}_{i}.\end{equation} \tag{ 30 }$$

Then they Taylor expand λ[σ] around σ = S.

$$\begin{equation}\frac{\mathrm{d}V}{\mathrm{d}t}=V\lambda [S]+\sum\limits _{i}{v}_{i}\left({\left.\frac{\mathrm{d}\lambda }{\mathrm{d}\sigma }\right\vert }_{S}({\sigma }_{\mathrm{b},i}-S)+\frac{1}{2}{\left.\frac{{\mathrm{d}}^{2}\lambda }{\mathrm{d}{\sigma }^{2}}\right\vert }_{S}{({\sigma }_{\mathrm{b},i}-S)}^{2}+\cdots \enspace \right).\end{equation} \tag{ 31 }$$

Using the fact that odd power terms of a should vanish due to symmetric cell division in terms of volume and an ansatz that the variance of the distribution of σ_b scales as a², Lin et al derived equation (29) without any knowledge about the distribution ψ[σ_b].

Equation (29) breaks down as a approaches 1, but when a is exactly 1, there exists a closed formula for Λ_p due to the self-similarity in a population tree as illustrated in figure 3. Lin et al showed that

$$\begin{equation}\sum\limits _{i=0}^{\infty }\mathrm{exp}\left(-{{\Lambda}}_{p}[a=1]\sum\limits _{j=0}^{i}\tau [jS]\right)=1.\end{equation} \tag{ 32 }$$

Zoom In Zoom Out Reset image size

Interpolating the solutions at a ≈ 0 and a = 1 with the minimal number of extrema, Lin et al infers the value of a_c (i.e. what value of a maximizes Λ_p for a given S) as described in figure 4. They found that the location of a_c depends on the relative magnitude of two critical values of S: S_c and S_*. Each of these two critical values arises from the solutions at a ≈ 0 and a = 1. S_c is where the second derivative of λ[S] is zero, which means that the fitness near a = 0 changes its sign when S = S_c. S_* is where Λ_p [a = 1] = Λ_p [a = 0]. See figure 4 for the inferred shape of the fitness curve at various values of S relative to S_c and S_*. Note that a_c is between 0 and 1 (marked by star in figure 4) in the case of beneficial key protein when S is between S_* and S_c.

Zoom In Zoom Out Reset image size

Similar to the previous paper by Lin et al, we use an interpolation to infer where the phase transition of a_c is. In addition, we are equipped with a more rigorous way of finding the population growth rate near a = 0 as described in section 3.4. We also have a transcendental equation for the slope of Λ_p [a] at a = 1, which we discuss in depth in appendix B. We denote the solution to equation (32) as Λ_p1 and expand Λ_p [a] = Λ_p1(1 + C₁(1 − a) + O((1 − a)²)). Finding the sign of C₁ makes it possible to predict whether Λ_p [a = 1] is a local minimum (C₁ > 0) or a local maximum (C₁ < 0). Using a similar moment expansion approach as we do for a ≈ 0, we find the first coefficient C₁ in appendix B as

$$\begin{equation}{C}_{1}=\frac{S}{2}\frac{{\sum }_{i=0}^{\infty }{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=0}^{i}\tau [jS]}{\sum }_{k=0}^{i}{\tau }^{\prime }[kS]\left({\sum }_{m=0}^{\infty }(m+1){\mathrm{e}}^{-{{\Lambda}}_{p}{\sum }_{n=0}^{m}\tau [nS]}-\frac{k(k+1)}{2}\right)}{{\sum }_{i=0}^{\infty }{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=0}^{i}\tau [jS]}{\sum }_{k=0}^{i}\tau [kS]}.\end{equation} \tag{ 33 }$$

Here we can define another critical value ${S}_{\mathrm{c}}^{a=1}$ which is the value of S that makes C₁ = 0. Because τ'[kS] > 0  (τ'[kS] < 0) for the damage (benefit) case, C₁ starts out positive (negative) for a small S and becomes negative (positive) as S exceeds ${S}_{\mathrm{c}}^{a=1}$. In other words, when $S{< }{S}_{\mathrm{c}}^{a=1}$, a = 1 is a local minimum (maximum) of Λ_p [a], whereas it becomes a local maximum (minimum) when $S{ >}{S}_{\mathrm{c}}^{a=1}$ in the damage (benefit) case. This provides us an additional piece of information that we can use to infer whether a_c = 1 for a given S. For instance, in the case of damage (λ₀ = 1, λ₁ = 0), ${S}_{\mathrm{c}}^{a=1}\approx 0.592$ is greater than S_* ≈ 0.558. This means that if S is between S_* and ${S}_{\mathrm{c}}^{a=1}$, Λ_p [a] should have a local maximum below a = 1. When S is above both S_* and ${S}_{\mathrm{c}}^{a=1}$, the population growth rate is maximized at a = 1 (i.e. a_c = 1). This means that the transition of a_c is not a single discrete jump from a_c = 0 to a_c = 1, but rather a jump to an intermediate a_c ∈ (0, 1) and then a smooth increase to a_c = 1. This process is illustrated in figure 5.

Zoom In Zoom Out Reset image size

In the case of beneficial key proteins, Lin et al found that a_c = 1 for small S and decreases smoothly to a_c = 0 as S increases, but the start of that transition was not specified. Based on appendix B, we find that a_c starts deviating from one at ${S}_{\mathrm{c}}^{a=1}=0.214$ (blue dotted line in figure 6(a)), consistent with the numerical solution of equation (10). Comparing the solid blue line (a_c predicted by the moment expansion to the sixth order) and the orange dots (numerical simulation), note that the moment expansion method does not successfully predict a_c when the true a_c value is close to 1. However, as S increases, and a_c becomes closer to 0, the answers agree better, since a² ≈ 0 our approximation holds.

Zoom In Zoom Out Reset image size

The moment expansion method can be easily generalized to the stochastic model where a is drawn from a distribution g(a) every time a cell divides. Equation (18) becomes

$$\begin{equation}\langle {\Delta}{\sigma }_{\mathrm{b}}^{n}\rangle =\frac{1}{{2}^{n-1}}\sum\limits _{m=0}^{\lfloor n/2\rfloor }\left(\genfrac{}{}{0pt}{}{n}{2m}\right)\int g(a){a}^{2m}\enspace \mathrm{d}a\left\langle {\Delta}{\sigma }_{\mathrm{b}}^{n-2m}{({\Delta}{\sigma }_{\mathrm{b}}+2S)}^{2m}\enspace {\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\Delta}{\sigma }_{\mathrm{b}}+S]}\right\rangle .\end{equation} \tag{ 34 }$$

This implies that in the small a limit, we can define an effective a, ${a}_{\text{eff}}=\sqrt{\langle {a}^{2}\rangle }$, for which the population growth rate Λ_p has the same formula independent of g(a).

$$\begin{equation}{{\Lambda}}_{p}=\lambda [S]\left(1+\frac{{\lambda }^{{\prime\prime}}[S]}{4\enspace \mathrm{ln}\enspace 2\lambda [S]}{S}^{2}{a}_{\text{eff}}^{2}\right).\end{equation} \tag{ 35 }$$

In figure 7, we numerically solve equation (10) where g(a) is a uniform distribution ($g(a)=U(\overline{a}-{\sigma }_{a},\overline{a}+{\sigma }_{\mathrm{b}})$ with $\overline{a}=0$ and σ_a = 0, 0.1, 0.2, ..., 1) and plot the population growth rate against ${a}_{\text{eff}}=\sqrt{\left\langle {a}^{2}\right\rangle }$. As predicted, a_eff behaves exactly as a does in the deterministic model for a broad range of σ_a . For a relatively large S and large σ_a , the population growth rate obtained by simulating the exponentially growing population starts to deviate from the second order prediction as shown in figure 7. This conclusion that stochastic segregation of key proteins increases effective asymmetry is analogous to the finding made by Chao et al [11] about damage partitioning in E. coli that larger fluctuations result in similar increase of the population growth as a fixed asymmetric ratio.

Zoom In Zoom Out Reset image size

So far, we have considered a model where the segregation ratio of key proteins is independent from the amount of key proteins each cell has. However, there is experimental evidence that the degree of asymmetry may depend on how much key proteins a cell has [2, 7]. Vedel et al [5] studied how protein aggregates in E. coli created in responds to heatshock are segregated differently between two daughters. One of their results is that the asymmetric segregation parameter a is a function of the amount of protein aggregates D, and if a is assumed to be constant, the fitness is underestimated. We can test the effect of dependence of a to D by replacing a with a[D] in equation (10).

$$\begin{equation}{\psi }_{\text{birth}}[{\sigma }_{\mathrm{b}}]=2\int {\psi }_{\text{birth}}[{\sigma }_{\mathrm{b}}^{\prime }]f({\sigma }_{\mathrm{b}}\vert {\sigma }_{\mathrm{b}}^{\prime },a[{\sigma }_{\mathrm{b}}^{\prime }+S]){\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\sigma }_{\mathrm{b}}^{\prime }]}\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }.\end{equation} \tag{ 36 }$$

Here we replaced D with σ_b' + S since we assume the cell volume always is 1 at birth and 2 at division. In figure 8, we show the numerical value of population growth rate increases due to asymmetry (Λ_p [a] − Λ_p [a = 0]), while a is a function of the amount of protein at division. Specifically, we set a[D] = tanh(kD) and tune k. In other words, we assume a to increase from zero initially with the slope k and asymptotically approach a = 1. Qualitatively, we expect the effective a to be small when S is small and to become closer to 1 as S gets larger. We can check this intuition by writing down the equation for the moments of ψ_birth analogous to equation (18).

$$\begin{equation}\langle {\Delta}{\sigma }_{\mathrm{b}}^{n}\rangle =\frac{1}{{2}^{n-1}}\sum\limits _{m=0}^{\lfloor n/2\rfloor }\left(\genfrac{}{}{0pt}{}{n}{2m}\right)\left\langle a{[{\Delta}{\sigma }_{\mathrm{b}}+2S]}^{2m}{\Delta}{\sigma }_{\mathrm{b}}^{n-2m}{({\Delta}{\sigma }_{\mathrm{b}}+2S)}^{2m}\enspace {\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\Delta}{\sigma }_{\mathrm{b}}+S]}\right\rangle .\end{equation} \tag{ 37 }$$

Zoom In Zoom Out Reset image size

If the slope of a[D] is small compared to a², we can approximate $a{[{\Delta}{\sigma }_{\mathrm{b}}+2S]}^{2m}\left.\right]\approx a[2S]$. Note that the population growth rate can be higher with an adaptive a (e.g. with k = 0.9, 0.55 < S < 0.65) than when a is fixed at a_c[S] (orange dots in figure 8). Qualitatively, this suggests that there is an advantage if each cell tunes the asymmetry of segregation according to its level of damage.

In this paper, we investigate how asymmetric segregation of key proteins in microbial organisms impact the population growth rate with a self-consistent equation (equation (10)) solved analytically and numerically. We adapt a physiological model proposed by Lin et al [8] and apply the transport equation to set up an equation for the distribution of protein concentration at birth ψ[σ_b] and population growth rate Λ_p (equation (10)). The equation can be solved by updating the solutions iteratively as elaborated in section 3.3. When the asymmetry is small (a ≪ 1), Λ_p can be approximated to a²ⁿ precision by solving a set of linear equations with respect to the moments of ψ[σ_b] and Λ_p . Similar approximations can be made from the opposite limit, a ≈ 1 (see appendix B). Changing the accumulation rate of key proteins, S, we confirm the sharp phase transition of a_c for a deleterious protein and a smooth transition for a beneficial protein as Lin et al reported [8]. In addition, since numerically solving equation (10) provides more precise prediction of Λ_p than a direct simulation of an exponentially growing population, we are able to find that the phase transition in the case of damage is composed of a discrete jump to a_c between 0 and 1 and a smooth increase to a_c = 1. We also explore different physiological models motivated by experimental studies. We discuss the case where the cell tunes the asymmetry with response to the amount of key proteins it has at division and how the adaptive strategy may minimize the disadvantage of asymmetry when the accumulation rate of damage is low. Finally, we show by generalizing our model that even if the segregation ratio is on average symmetric, the noise can have a similar effect as an asymmetric segregation.

There are many other aspects of asymmetric division in biology where our theoretical framework might be useful for building a theoretical model [1, 16]. First, there are many other ways that microbes develop phenotypic diversity rather than due to asymmetric key protein segregation. For instance, budding yeast (Saccharomyces cerevisiae) has cell size difference between a daughter (smaller cell) and mother cell (bigger cell), whereas when a Caulobacter crescentus cell divides, one is non-motile and reproductive and the other swarms and does not divide until maturation [17, 18]. In the former case, introduction of volume asymmetry turns out to create interesting correlations between the noise in growth and the population growth rate [19]. Second, mortality rate rather than volumetric growth rate may depend on the key protein concentration [20]. Lastly, connecting back to multicellular organisms, where the theory of aging (soma theory) was first proposed [21], can be interesting both mathematically and biologically. Major difference from microbial system is that the cells interact with the neighbors and use the cues to make decisions [22]. Considering the geometry of the tissue might complicate the model beyond the level of a self-consistent equation can capture, but plants which have rigid tissue and simpler geometry may be a good starting point for a theoretical investigation [23].

The authors thank the Amir group, the Desai lab, Andrew Murray, Jane Kondev and Doeke Hekstra for helpful discussions. J M acknowledges Yubo Su for helping improve numerical simulations. A A was supported by the NSF CAREER Award Number 1752024.

When we explicitly simulate the exponentially growing population, we start with 100 cells with uniform distribution of age. We run until the number of cells reaches 5 × 10⁶ and find the slope of ln(∑V(t)) after t = 5. This use of sum of volume rather than the number of cells for finding the population growth rate is because the volume increases continuously and the growth rate should be equivalent if the cell size is controlled. To solve equation (10), we discretize the range of σ_b and repeat the process described in section 3.3 until the difference between 1 and $2\int \int {\psi }_{\text{birth}}^{(n)}[{\sigma }_{\mathrm{b}}^{\prime }]g({a}^{\prime }){\mathrm{e}}^{-{{\Lambda}}_{p}^{(n-1)}\tau [{\sigma }_{\mathrm{b}}^{\prime }]}\mathrm{d}{a}^{\prime }\enspace \mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }$ is smaller than some tolerance value. This implies that ${{\Lambda}}_{p}^{(n-1)}\approx {{\Lambda}}_{p}^{(n)}$ and our solution for Λ_p has converged. Python scripts for the simulations are on our Github (https://github.com/jiseonmin/asymmetric_segregation).

The challenge for a highly asymmetric case is that the distribution ψ_birth[σ_b] is not symmetric around S. Because one of two newborn cells is rejuvenated at each division, we would expect the distribution to be skewed to the right and maximized near σ_b = 0. This means that we should expand ψ_birth[σ_b] around different values depending on when the cell has been rejuvenated for the last time (i.e. generational age), in order to solve the self-consistent equation. Categorizing cells based on their generational age, we can rewrite the self-consistent equation (10) as,

$$\begin{align}{\psi }_{0}[{\sigma }_{\mathrm{b}}]& =\int \sum\limits _{i=0}^{\infty }{\psi }_{i}[{\sigma }_{\mathrm{b}}^{\prime }]f({\sigma }_{\mathrm{b}},\mathrm{R}\vert {\sigma }_{\mathrm{b}}^{\prime }){\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\sigma }_{\mathrm{b}}^{\prime }]}\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime },\\ {\psi }_{i}[{\sigma }_{\mathrm{b}}]& =\int {\psi }_{i-1}[{\sigma }_{\mathrm{b}}^{\prime }]f({\sigma }_{\mathrm{b}},\mathrm{S}\vert {\sigma }_{\mathrm{b}}^{\prime }){\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\sigma }_{\mathrm{b}}^{\prime }]}\enspace \mathrm{d}{\sigma }_{\mathrm{b}}\quad (i{ >}0).\end{align} \tag{ B.1 }$$

Here, ψ_i [σ_b] is the probability that a cell has protein concentration σ_b at birth and the generational age i. The transition function for the rejuvenated cell is $f(\sigma ,\mathrm{R}\vert {\sigma }_{\mathrm{b}}^{\prime })=\delta \left(\sigma -\frac{1-a}{2}({\sigma }_{\mathrm{b}}^{\prime }+S)\right)$, and the senescent is $f(\sigma ,\mathrm{S}\vert {\sigma }_{\mathrm{b}}^{\prime })=\delta \left(\sigma -\frac{1+a}{2}({\sigma }_{\mathrm{b}}^{\prime }+S)\right)$. Note that the rejuvenating transition function does not depend on σ_b' (parent cell) if a = 1. This is equivalent to the property used by Lin et al [8] to find a closed formula for Λ_p [a = 1]: there is a self-symmetry in a population tree in a stable exponential growth.

$$\begin{equation}1=\sum\limits _{i=0}^{\infty }\mathrm{exp}\left[-\sum\limits _{j=0}^{i}{{\Lambda}}_{p}\tau [jS]\right].\end{equation} \tag{ B.2 }$$

Here, we can use the fact that every ψ_i [σ_b] is a delta function if a = 1.

$$\begin{align}{\psi }_{i}[{\sigma }_{\mathrm{b}}]& =\int {\psi }_{i-1}[{\sigma }_{\mathrm{b}}^{\prime }]\delta ({\sigma }_{\mathrm{b}}-({\sigma }_{\mathrm{b}}^{\prime }+S)){\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\sigma }_{\mathrm{b}}^{\prime }]}\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }\\ & ={\psi }_{i-1}[{\sigma }_{\mathrm{b}}-S]{\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\sigma }_{\mathrm{b}}-S]}\\ & ={\psi }_{0}[{\sigma }_{\mathrm{b}}-iS]{\mathrm{e}}^{-{\sum }_{j=1}^{i}{{\Lambda}}_{p}\tau [{\sigma }_{\mathrm{b}}-jS]}\quad (i{ >}0),\end{align} \tag{ B.3 }$$

$$\begin{align}{\psi }_{0}[{\sigma }_{\mathrm{b}}]& =\int \sum\limits _{i=0}^{\infty }{\psi }_{i}[{\sigma }_{\mathrm{b}}^{\prime }]\delta ({\sigma }_{\mathrm{b}}){\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\sigma }_{\mathrm{b}}^{\prime }]}\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }\\ & =\delta ({\sigma }_{\mathrm{b}})\int \sum\limits _{i=0}^{\infty }{\psi }_{0}[{\sigma }_{\mathrm{b}}-iS]\mathrm{exp}\left[-\sum\limits _{j=0}^{i}{{\Lambda}}_{p}\tau [{\sigma }_{\mathrm{b}}-jS]\right]\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }.\end{align} \tag{ B.4 }$$

By integrating ψ₀[σ_b] with respect to σ_b, we get

$$\begin{equation}1=\sum\limits _{i=0}^{\infty }\mathrm{exp}\left[-\sum\limits _{j=0}^{i}{{\Lambda}}_{p}\tau \enspace [(i-j)S]\right]=\sum\limits _{i=0}^{\infty }\mathrm{exp}\left[-\sum\limits _{j=0}^{i}{{\Lambda}}_{p}\tau [jS]\right],\end{equation} \tag{ B.5 }$$

which is consistent with the previous paper [8]. Let us call the unique solution for equation (B.5), Λ_p,1. Near a = 1, we can find the approximate population growth rate by expanding Λ_p = Λ_p1(1 + C₁(1 − a) + C₂(1 − a)² + ⋯) where a < 1. As a deviates from 1, we expect each delta function to broaden near σ_b = iS. Similar to the case of a small a, we will find the moments of ψ_i [σ_b] starting from the lowest order. Here, we will investigate the first order coefficient C₁. Importantly, finding the sign of C₁ will allow us to predict whether Λ_p [a = 1] is a local minimum (C₁ > 0) or a local maximum (C₁ < 0).

Let us define δa = 1 − a. To find C₁ we need to solve equation (B.1) to O(δa). Thus, we can rewrite the equations as

$$\begin{align}{\psi }_{0}[{\sigma }_{\mathrm{b}}]& =\int \sum\limits _{i=0}^{\infty }{\psi }_{i}[{\sigma }_{\mathrm{b}}^{\prime }]\delta \left({\sigma }_{\mathrm{b}}-\frac{\delta a}{2}({\sigma }_{\mathrm{b}}^{\prime }+S)\right)\\ & \quad {\times}\mathrm{exp}\left[-{{\Lambda}}_{p1}\tau [iS](1+{C}_{1}\delta a)\left(1+\frac{{\tau }^{\prime }[iS]}{\tau [iS]}({\sigma }_{\mathrm{b}}^{\prime }-iS)\right)\right]\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime },\\ {\psi }_{i}[{\sigma }_{\mathrm{b}}]& =\int {\psi }_{i-1}[{\sigma }_{\mathrm{b}}^{\prime }]\delta \left({\sigma }_{\mathrm{b}}-(1-\delta a/2)({\sigma }_{\mathrm{b}}^{\prime }+S)\right)\\ & \quad {\times}\mathrm{exp}\left[-{{\Lambda}}_{p1}\tau [(i-1)S](1+{C}_{1}\delta a)\left(1+\frac{{\tau }^{\prime }[(i-1)S]}{\tau [(i-1)S]}({\sigma }_{\mathrm{b}}^{\prime }-(i-1)S)\right)\right]\mathrm{d}{\sigma }_{\mathrm{b}}.\end{align} \tag{ B.6 }$$

For brevity, from now on, we shall write Ψ_i = ∫ψ_i [σ_b]dσ_b and ${\left\langle F[\delta {\sigma }_{\mathrm{b}}]\right\rangle }_{i}=\int {\psi }_{i}[{\sigma }_{\mathrm{b}}]F[{\sigma }_{\mathrm{b}}-iS]\mathrm{d}{\sigma }_{\mathrm{b}}/{{\Psi}}_{i}$. First, let us integrate the ith equation (i > 0),

$$\begin{align}{{\Psi}}_{i}& =\int {\psi }_{i-1}[{\sigma }_{\mathrm{b}}^{\prime }]{\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\sigma }_{\mathrm{b}}^{\prime }]}\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }\\ & ={\mathrm{e}}^{-{{\Lambda}}_{p1}\tau [(i-1)S]}\int {\psi }_{i-1}[{\sigma }_{\mathrm{b}}^{\prime }]\left(1-{{\Lambda}}_{p1}\tau [(i-1)S]{C}_{1}\delta a-{{\Lambda}}_{p1}{\tau }^{\prime }[(i-1)S]({\sigma }_{\mathrm{b}}^{\prime }-(i-1)S)\right)\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }\\ & ={{\Psi}}_{i-1}{\mathrm{e}}^{-{{\Lambda}}_{p1}\tau [(i-1)S]}\left(1-{{\Lambda}}_{p1}\tau [(i-1)S]{C}_{1}\delta a-{\left\langle \delta {\sigma }_{\mathrm{b}}\right\rangle }_{i-1}{{\Lambda}}_{p1}{\tau }^{\prime }[(i-1)S]\right)\\ & ={{\Psi}}_{0}{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=0}^{i-1}\tau [jS]}\left(1-{{\Lambda}}_{p1}\sum\limits _{k=0}^{i-1}\tau [kS]{C}_{1}\delta a-{\left\langle \delta {\sigma }_{\mathrm{b}}\right\rangle }_{k}{{\Lambda}}_{p1}{\tau }^{\prime }[kS]\right)+O(\delta {a}^{2}).\end{align} \tag{ B.7 }$$

For i = 0 we have

$$\begin{align}& {{\Psi}}_{0}=\sum\limits _{i=0}^{\infty }\int {\psi }_{i}[{\sigma }_{\mathrm{b}}^{\prime }]{\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\sigma }_{\mathrm{b}}^{\prime }]}\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }\\ & \quad =\sum\limits _{i=1}^{\infty }{{\Psi}}_{i}\\ & \quad ={{\Psi}}_{0}\sum\limits _{i=0}^{\infty }{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=0}^{i}\tau [jS]}\left(1-{{\Lambda}}_{p1}\sum\limits _{k=0}^{i}\tau [kS]{C}_{1}\delta a-{\left\langle \delta {\sigma }_{\mathrm{b}}\right\rangle }_{k}{{\Lambda}}_{p1}{\tau }^{\prime }[kS]\right),\end{align} \tag{ B.8 }$$

to O(δa) accuracy. Using the definition of Λ_p1 − $1={\sum }_{i=0}^{\infty }{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=0}^{i}\tau [jS]}$, we can extract an expression for the coefficient C₁.

$$\begin{align}{C}_{1}& =-\frac{{\sum }_{i=0}^{\infty }{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=0}^{i}\tau [jS]}{\sum }_{k=0}^{i}{\tau }^{\prime }[kS]{\langle \delta {\sigma }_{\mathrm{b}}\rangle }_{k}}{{\sum }_{i=0}^{\infty }{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=0}^{i}\tau [jS]}{\sum }_{k=0}^{i}\tau [kS]\delta a}\\ & =-\frac{{\sum }_{k=0}^{\infty }{\tau }^{\prime }[kS]{\left\langle \delta {\sigma }_{\mathrm{b}}\right\rangle }_{k}{\sum }_{i=k}^{\infty }{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=0}^{i}\tau [jS]}}{\delta a{\sum }_{k=0}^{\infty }\tau [kS]{\sum }_{i=k}^{\infty }{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=0}^{i}\tau [jS]}}.\end{align} \tag{ B.9 }$$

Let us find ${\left\langle \delta {\sigma }_{\mathrm{b}}\right\rangle }_{i}$ for i > 0.

$$\begin{align}{\left\langle {\sigma }_{\mathrm{b}}\right\rangle }_{i}& =\frac{1}{{{\Psi}}_{i}}{\mathrm{e}}^{-{{\Lambda}}_{p1}\tau [(i-1)S]}\int {\psi }_{i-1}[{\sigma }_{\mathrm{b}}^{\prime }]\left(\left(1-\frac{\delta a}{2}\right)\left({\sigma }_{\mathrm{b}}^{\prime }-(i-1)S+iS\right)-iS\right)\\ & \quad {\times}\left(1-{{\Lambda}}_{p1}\tau [(i-1)S]{C}_{1}\delta a-{{\Lambda}}_{p1}{\tau }^{\prime }[(i-1)S]({\sigma }_{\mathrm{b}}^{\prime }-(i-1)S)\right)\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }\\ & =-\frac{\delta a}{2}iS+\frac{1}{{{\Psi}}_{i}}{\mathrm{e}}^{-{{\Lambda}}_{p1}\tau [(i-1)S]}{{\Psi}}_{i-1}{\left\langle \delta {\sigma }_{\mathrm{b}}\right\rangle }_{i-1}\\ & =-\frac{\delta a}{2}\sum\limits _{k=0}^{i}(i-k)S\frac{{{\Psi}}_{i-k}}{{{\Psi}}_{i}}{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=1}^{k}\tau [(i-j)]}+\frac{{{\Psi}}_{0}}{{{\Psi}}_{i}}{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=0}^{i-1}\tau [jS]}{\left\langle \delta {\sigma }_{\mathrm{b}}\right\rangle }_{0}\\ & =-\delta a\frac{S}{2}\frac{i(i+1)}{2}+{\left\langle \delta {\sigma }_{\mathrm{b}}\right\rangle }_{0}+O(\delta {a}^{2}).\end{align} \tag{ B.10 }$$

If i = 0,

$$\begin{align}{\left\langle \delta {\sigma }_{\mathrm{b}}\right\rangle }_{0}& =\frac{1}{{{\Psi}}_{0}}\frac{\delta a}{2}\sum\limits _{m=0}^{\infty }{\psi }_{m}[{\sigma }_{\mathrm{b}}^{\prime }]\left({\sigma }_{\mathrm{b}}^{\prime }-mS+(m+1)S\right)\\ & \quad {\times}\left[1-{{\Lambda}}_{p1}\tau [mS]{C}_{1}\delta a-{{\Lambda}}_{p1}{\tau }^{\prime }[mS]({\sigma }_{\mathrm{b}}^{\prime }-mS)\right]\mathrm{d}{\sigma }_{\mathrm{b}}^{\prime }\\ & =\frac{1}{{{\Psi}}_{0}}\frac{\delta a}{2}\sum\limits _{m=0}^{\infty }(m+1)S{{\Psi}}_{m+1}\end{align} \tag{ B.11 }$$

$$\begin{equation}\qquad =\frac{\delta a}{2}\sum\limits _{m=0}^{\infty }(m+1)S{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{n=0}^{m}\tau [nS]}+O(\delta {a}^{2}).\end{equation} \tag{ B.12 }$$

This shows that ${\left\langle \delta {\sigma }_{\mathrm{b}}\right\rangle }_{i}$ is positive for small i and becomes more and more negative as i increases. Applying this to equation (B.9), we can find the slope near a = 1.

$$\begin{equation}{C}_{1}=\frac{S}{2}\frac{{\sum }_{i=0}^{\infty }{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=0}^{i}\tau [jS]}{\sum }_{k=0}^{i}{\tau }^{\prime }[kS]\left(-\frac{k(k+1)}{2}+{\sum }_{m=0}^{\infty }(m+1){\mathrm{e}}^{-{{\Lambda}}_{p}{\sum }_{n=0}^{m}\tau [nS]}\right)}{{\sum }_{i=0}^{\infty }{\mathrm{e}}^{-{{\Lambda}}_{p1}{\sum }_{j=0}^{i}\tau [jS]}{\sum }_{k=0}^{i}\tau [kS]}.\end{equation} \tag{ B.13 }$$

In order to calculate Λ_p [a] to O(a²ⁿ ) accuracy, we write Λ_p and τ[σ_b] as

$$\begin{equation}{{\Lambda}}_{p}=\frac{\mathrm{ln}\enspace 2}{\tau [S]}\left(1+\sum\limits _{k=1}^{n}{C}_{2k}{a}^{2k}\right),\end{equation} \tag{ C.1 }$$

$$\begin{equation}\tau [{\sigma }_{\mathrm{b}}^{\prime }]=\tau [S]\sum\limits _{i=0}^{2n}\frac{{\tau }^{(i)}[S]}{i!\tau [S]}{\Delta}{\sigma }_{\mathrm{b}}^{\prime i}.\end{equation} \tag{ C.2 }$$

Using these, we write

$$\begin{equation}{\mathrm{e}}^{-{{\Lambda}}_{p}\tau [{\sigma }_{\mathrm{b}}]}=\mathrm{exp}\left[-\mathrm{ln}\enspace 2\left(1+\sum\limits _{k=1}^{n}{C}_{2k}{a}^{2k}\right)\sum\limits _{i=0}^{2n}\frac{{\tau }^{(i)}[S]}{i!\tau [S]}{\Delta}{\sigma }_{\mathrm{b}}^{i}\right].\end{equation} \tag{ C.3 }$$

Next, we should expand equation (18) and truncate any terms with higher order than a²ⁿ . We can also rewrite $\langle {\sigma }_{\mathrm{b}}^{m}\rangle \enspace (1{\leqslant}m{\leqslant}2n)$ as

$$\begin{equation}\langle {\Delta}{\sigma }_{\mathrm{b}}^{2k}\rangle =\sum\limits _{j=k}^{n}{D}_{2k}^{(2j)}{a}^{2j},\end{equation} \tag{ C.4 }$$

$$\begin{equation}\langle {\Delta}{\sigma }_{\mathrm{b}}^{2k-1}\rangle =\sum\limits _{j=k}^{n}{D}_{2k+1}^{(2j)}{a}^{2j}.\end{equation} \tag{ C.5 }$$

Maintaining only terms to O(a²ⁿ ) and matching the coefficients of a^2j (1 ⩽ j ⩽ n), we can solve the system of equations for D and C's. We provide the SymPy code to find the coefficients on our Github repository (https://github.com/jiseonmin/asymmetric_segregation).

Here we present some examples of the distribution of protein concentration at birth ψ_birth[σ_b], which we find by numerically solving equation (10) (red circles). In each simulation, we compare the numerically found distribution to a log-normal distribution and a normal distribution that have the same mean and variance as the numerical solution. In other words, we test whether it is possible to approximate the distribution of key protein concentration to a log-normal or a normal distribution and simplify equation (10) (figures D1–D3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When a is small and the segregation is noisy (σ_a > 0), the distribution can be approximated as a Gaussian distribution or a log-normal distribution. However, as a increases, the skewness of the distribution grows and the distribution becomes multi-modal. This implies that with small average asymmetry, the first two moments of the distribution have sufficient information for approximately solving equation (10) for Λ_p , but as the segregation becomes more and more asymmetric, higher order terms of a are necessary to accurately calculate the population growth rate. The exact value of a where the Gaussian approximation breaks down depends on the shape of λ[σ] and S.