We propose a stochastic model of cell division, death and differentiation across an ordered sequence of compartments. Cells in a given compartment may represent a common spatial location within the body, or a common phenotype (i.e., representing a biological state, defined by morphology and/or function). In practice, this means that cells in a compartment behave equally, with the same per-cell rates for a given division, death, or differentiation event. We consider a sequence of compartments \(C_i\), \(i\in \{1,\ldots ,N\}\), which cells follow, behaving independently of each other, as generally assumed in the theory of continuous-time branching processes40.Our general stochastic model considers a number of cellular events inspired by some recent mathematical models3,9,14,42. We assume that each of these events takes place at a given per-cell rate. Particular situations of interest arise from setting some of the rates equal to zero, so that those events cannot take place, as we illustrate for three different case studies in Section “Results”. Each cell in a given compartment, \(C_i\), can divide, die or exit to one of the two adjacent compartments. When a division occurs, both daughter cells might belong to the same compartment as the mother (this event is referred to as self-renewal), both daughter cells might instantaneously move to the next compartment (symmetric division), or one daughter cell might belong to the same compartment as the mother, and the other to the next compartment (asymmetric division)3,9.Figure 1General stochastic model of cell division, death and differentiation for an ordered sequence of compartments. Grey cells represent a death event. Self-renewal events take place with rate \(\lambda_i\), symmetric division events with rate \(s_i\), and asymmetric division events with rate \(a_i\). Differentiation events happen with per-cell rate \(\nu_i\) (forward) or \(\xi_i\) (backward). Death events have per-cell rate \(\mu_i\). Each per-cell rate is a real positive number.The stochastic model, shown in Fig. 1, is a continuous-time Markov chain (CTMC) \({{{\mathcal {X}}}}=\{({\varvec{C}}_1(t),{\varvec{C}}_2(t),\ldots ,{\varvec{C}}_N(t)):\ t\ge 0\}\), where \({\varvec{C}}_i(t)\) is a random variable that represents the number of cells in compartment \(C_i\) at time t, with state space given by \(\mathcal{S}=\{0,1,2,\ldots \}^N={\mathbb {N}}_0^N\). At any given time, the process can be at a particular state \((n_1,\ldots ,n_N)\in {{{\mathcal {S}}}}\), so that there are \(n_i\) cells in compartment \(C_i\). Thus, cell events, labelled E1 to E5 below, represent transitions between states as follows:
(E1)
Self-renewal (cell division where both daughter cells remain in the same compartment as the mother) can occur in any compartment \(C_i\), with per-cell rate \(\lambda _i\), for \(i\in \{1,\ldots ,N\}\), $$\begin{aligned} (n_1,\ldots ,n_{i-1},n_i,n_{i+1},\ldots ,n_N)\rightarrow_{\lambda_i n_i}\rightarrow & ~ (n_1,\ldots ,n_{i-1},n_i+1,n_{i+1},\ldots ,n_N) \nonumber . \end{aligned}$$ The previous notation indicates that the event moving the process from state \((n_1,\ldots ,n_{i-1},n_i,n_{i+1},\ldots ,n_N)\) to state \((n_1,\ldots ,n_{i-1},n_i+1,n_{i+1},\ldots ,n_N)\) happens at rate \(\lambda _i n_i\), and represents self-renewal.
(E2)
Symmetric division (cell division where both daughter cells instantaneously move to the next compartment) can occur in compartment \(C_i\) with per-cell rate \(s_i\), for \(i\in \{1,\ldots ,N-1\}\), $$(n_{1} , \ldots ,n_{{i – 1}} ,n_{i} ,n_{{i + 1}} , \ldots ,n_{N} ) \rightarrow_{s_i n_i}\rightarrow ~ (n_{1} , \ldots ,n_{{i – 1}} ,n_{i} – 1,n_{{i + 1}} + 2, \ldots ,n_{N} ).$$
(E3)
Asymmetric division (cell division where one of the daughter cells remains in the same compartment as the mother, while the other goes to the next compartment) can occur in compartment \(C_i\) with per-cell rate \(a_i\), for \(i\in \{1,\ldots ,N-1\}\), $$\begin{aligned} (n_1,\ldots ,n_{i-1},n_i,n_{i+1},\ldots ,n_N) \rightarrow_{a_i n_i}\rightarrow & ~ (n_1,\ldots ,n_{i-1},n_i,n_{i+1}+1,\ldots ,n_N). \end{aligned}$$
(E4)
Differentiation (or migration) between adjacent compartments can occur with per-cell rate \(\nu _i\), for \(i\in \{1,\ldots ,N-1\}\) and \(\xi _i\), for \(i\in \{2,\ldots ,N\}\), $$\begin{aligned} (n_1,\ldots ,n_{i-1},n_i,n_{i+1},\ldots ,n_N) \rightarrow_{\nu_i n_i}\rightarrow & ~ (n_1,\ldots ,n_{i-1},n_i-1,n_{i+1}+1,\ldots ,n_N),\\ (n_1,\ldots ,n_{i-1},n_i,n_{i+1},\ldots ,n_N) \rightarrow_{\xi_i n_i}\rightarrow & ~ (n_1,\ldots ,n_{i-1}+1,n_i-1,n_{i+1},\ldots ,n_N). \end{aligned}$$
(E5)
Cells can die in any compartment \(C_i\) with per-cell rate \(\mu _i\), \(i\in \{1,\ldots ,N\}\), $$\begin{aligned} (n_1,\ldots ,n_{i-1},n_i,n_{i+1},\ldots ,n_N) \rightarrow_{\mu_i n_i}\rightarrow & ~ (n_1,\ldots ,n_{i-1},n_i-1,n_{i+1},\ldots ,n_N). \end{aligned}$$
We assume that cells in the last compartment, \(C_N\), cannot symmetrically or asymmetrically divide, or differentiate to the next compartment. All rates are assumed to be positive real numbers.Mean number of cells in each compartmentWe first study the dynamics of the mean number of cells in each compartment at time t, \({\mathbb {E}}[{\varvec{C}}_{i}(t)]\), described by the following system of ordinary differential equations44$$\begin{aligned} \frac{d\, {\mathbb {E}}[{\varvec{C}}_{1}(t)]}{dt}= & {} -(\mu _1+\nu _1+s_1-\lambda _1) {\mathbb {E}}[{\varvec{C}}_{1}(t)] + \xi _2{\mathbb {E}}[{\varvec{C}}_{2}(t)], \nonumber \\ \frac{d\, {\mathbb {E}}[{\varvec{C}}_{i}(t)]}{dt}= & {} (\nu _{i-1}+a_{i-1}+2s_{i-1})\,{\mathbb {E}}[{\varvec{C}}_{i-1}(t)] – (\mu _i+\nu _i+\xi _i+s_i-\lambda _i) \, {\mathbb {E}}[{\varvec{C}}_{i}(t)]+\xi _{i+1}\, {\mathbb {E}}[{\varvec{C}}_{i+1}(t)], \quad i \in \{2,\ldots ,N-1\},\nonumber \\ \frac{d\, {\mathbb {E}}[{\varvec{C}}_{N}(t)]}{dt}= & {} (\nu _{N-1}+a_{N-1}+2s_{N-1})\, {\mathbb {E}}[{\varvec{C}}_{N-1}(t)]- (\mu _N+\xi _N-\lambda _N) {\mathbb {E}}[{\varvec{C}}_{N}(t)], \end{aligned}$$
(1)
where \({\mathbb {E}}[{\varvec{C}}_i(t)]\) represents the expectation of the random variable \({\varvec{C}}_i(t)\). These equations constitute a homogeneous first-order linear system of ODEs with constant coefficients, which can be written more succinctly in matrix form as follows$$\begin{aligned} \frac{d\textbf{C}(t)}{dt} = \textbf{A} \> \textbf{C}(t), \end{aligned}$$
(2)
where$$\begin{aligned} \textbf{C}(t) = \left( \begin{array}{c} {\mathbb {E}}[{\varvec{C}}_{1}(t)]\\ {\mathbb {E}}[{\varvec{C}}_{2}(t)]\\ \vdots \\ {\mathbb {E}}[{\varvec{C}}_{N-1}(t)]\\ {\mathbb {E}}[{\varvec{C}}_{N}(t)]\end{array}\right) ,\quad \textbf{A} = \left( \begin{array}{ccccc} -\Delta _1 &{} \xi _{2} &{} 0 &{} \cdots &{} 0\\ \Lambda _1 &{} -\Delta _2 &{} \xi _{3}&{} \ldots &{}\vdots \\ 0 &{}\ddots &{} \ddots &{}\ddots &{}0\\ \vdots &{} \vdots &{} \Lambda _{N-2} &{} -\Delta _{N-1}&{} \xi _{N}\\ 0 &{} \cdots &{} 0 &{} \Lambda _{N-1} &{}-\Delta _N \end{array} \right) , \end{aligned}$$
(3)
and$$\begin{aligned} \Delta _1= & {} \mu _1+\nu _1+ s_1 -\lambda _1, \\ \Delta _i= & {} \mu _i+\nu _i+s_i+ \xi _{i}-\lambda _i, \quad \Lambda _{i-1} \ =\ \nu _{i-1}+a_{i-1}+2s_{i-1},\ i\in \{2,\ldots ,N-1\}, \\ \Delta _N= & {} \mu _N+ \xi _{N} – \lambda _N,\quad \Lambda _{N-1} \ =\ \nu _{N-1}+a_{N-1}+2s_{N-1}. \end{aligned}$$The initial value problem of Eq. (2) with \(\textbf{C}_0 = \textbf{C}(0)\) has a unique solution [45, Theorem 3.9] given by$$\begin{aligned} \textbf{C}(t) = e^{\textbf{A}t}{} \textbf{C}_0, \end{aligned}$$
(4)
where \(e^{\textbf{A}t}\) represents the matrix exponential$$\begin{aligned} e^{\textbf{A}t} =\textbf{I}+\textbf{A}t+\textbf{A}^2 \frac{t^2}{2!}+\textbf{A}^3 \frac{t^3}{3!}+ \cdots = \sum _{i=0}^{+\infty }\frac{(\textbf{A}t)^i}{i!}. \end{aligned}$$The system of equations (2) admits \(\lim _{t\rightarrow + \infty }{} \textbf{C}(t)=\textbf{0}_N\) (column vector of zeros) as asymptotic solution. It is exponentially stable (solutions of the system for any initial condition converge to zero) if and only if each eigenvalue of \(\textbf{A}\) has a negative real part (see, for example, Ref.[45, Corollaries 3.5 and 3.6]).Figure 2Irreversible stochastic model of cell division, death and forward differentiation for an ordered sequence of compartments. Grey cells represent death events. Self-renewal events occur with per-cell rate \(\lambda_i\), symmetric division events with rate \(s_i\), and asymmetric division events with rate \(a_i\). Differentiation events happen with per-cell rate \(\nu_i\). Death events have (per-cell) rate \(\mu_i\). Each per-cell rate is a real positive number.We note that, for certain biological applications, some of the rates in Fig. 1 will be zero, and thus, the analysis of such systems would simplify. For instance, differentiation may be irreversible9,14, so that \(\xi _i=0\) for \(i\in \{2,\ldots ,N\}\) as shown in Fig. 2. We will refer to this scenario as the irreversible model. In this case, and if one considers a single progenitor cell starting in compartment \(C_1\) at time \(t=0\), \(\textbf{C}(0)=(1,0,\ldots ,0)^T\), it is possible to obtain the mean number of cells for any compartment. In fact, one has$$\begin{aligned} {\mathbb {E}}[{\varvec{C}}_{i}(t)]= & {} \left\{ \begin{array}{ll} e^{-\Delta _1t}, &{} i=1,\\ \left( \prod \limits _{l=1}^{i-1} \Lambda _l \right) \sum \limits _{j=1}^i e^{-\Delta _j t} \prod \limits _{\begin{array}{c} m=1 \\ m\ne j \end{array}}^i (\Delta _m-\Delta _j)^{-1},&i\in \{2,\ldots ,N\}. \end{array}\right. \end{aligned}$$
(5)
Equation (5) is well-defined if \(\Delta _i\ne \Delta _j\) for pairs (i, j) with \(i \ne j\). If this is not the case, alternative analytic solutions can be found. For example, if \(\Delta _i=\Delta _j\) for all \(i,j\in \{1,\ldots ,N\}\), Eq. (5) simplifies to$$\begin{aligned} {\mathbb {E}}[{\varvec{C}}_{i}(t)] = \left( \prod \limits _{l=1}^{i-1} \Lambda _l \right) \frac{t^{i-1}}{(i-1)!}e^{-\Delta _it},\quad t\ge 0. \end{aligned}$$
(6)
We note that this is consistent with existing results in the literature. In particular, the most general solution for the case of repeated eigenvalues can be found in Ref.46, which analysed the radioactive decay of a chain of nuclides. It is clear that in the irreversible model \(\lim _{t\rightarrow +\infty }{\mathbb {E}}[{\varvec{C}}_i(t)]=0\), when \(\Delta _i>0\) \(\forall i\in \{1,\ldots ,N\}\). This is consistent since \(\{-\Delta _i:\ i\in \{1,\ldots ,N\}\}\) are the eigenvalues of \(\textbf{A}\) in this case.For biological applications3,14 such as those considered in Section “Results”, it is of interest to determine the cumulative average number of cells that arrive to the final compartment, \(C_N\), starting with a single “progenitor” or precursor cell in compartment \(C_1\). To this end, in the irreversible model one can set \(\lambda _N=\mu _N=0\), so that cells arriving into \(C_N\) accumulate and can be counted. Then, for the last compartment \(\Delta _N=0\) and thus, Eq. (5) leads to$$\begin{aligned} {\mathbb {E}}[{\varvec{C}}_{N}(t)] = \left( \prod \limits _{l=1}^{N-1} \Lambda _l \right) \left[ \sum \limits _{j=1}^{N-1} e^{-\Delta _j t} \prod \limits _{\begin{array}{c} m=1 \\ m\ne j \end{array}}^{N-1} (\Delta _m-\Delta _j)^{-1} + \prod \limits _{\begin{array}{c} m=1 \end{array}}^{N-1} \Delta _m^{-1} \right] , \end{aligned}$$
(7)
which is well-defined if \(\Delta _i\ne \Delta _j\) for all \(i,j\in \{1,\ldots ,N-1\}\). From the previous equation, we have$$\begin{aligned} \lim _{t\rightarrow +\infty } {\mathbb {E}}[ {\varvec{C}}_{N}(t)]=\prod _{i=1}^{N-1}\frac{\Lambda _i}{\Delta _i}. \end{aligned}$$
(8)
Interestingly, this limit also holds if \(\Delta _i=\Delta _j\) for all \(i,j \in \{1,\ldots ,N-1\}\) and \(\Delta _N=0\). In this case one can write$$\begin{aligned} {\mathbb {E}}[{\varvec{C}}_{N}(t)]= & {} \prod _{i=1}^{N-1}\frac{\Lambda _i}{\Delta _i} -\sum _{j=1}^{N-1} {\mathbb {E}}[{\varvec{C}}_j(t)] \prod _{i=j}^{N-1}\frac{\Lambda _i}{\Delta _i}. \end{aligned}$$
(9)
Under population extinction conditions (that is, when \(\Delta _i>0\) for all \(i\in \{1,\ldots ,N-1\}\), so that the population of cells in intermediate compartments, \(\{1,\ldots ,N-1\}\), dies out and only “exiting” cells remain in \(C_N\) at late times), \(\lim _{t\rightarrow +\infty } {\mathbb {E}}[ {\varvec{C}}_{i}(t)]=0\) for \(i\in \{1,\ldots ,N-1\}\), and thus Eq. (8) also holds.A particular feature of this system is that cells behave independently from each other. This means that the dynamics of the progeny of a set of M progenitor (or precursor) cells in compartment \(C_1\) at time \(t=0\), can be analysed as M independent stochastic processes. Thus, in Section “The progeny of a single progenitor cell” we consider a number of summary statistics of interest related to the stochastic journeys of cell progenies from a single cell starting at a given compartment (typically compartment \(C_1\)).The progeny of a single progenitor cellFor a cell starting in compartment \(C_i\), we can define \(G_i\) to be the random variable representing the total number of cells in the progeny of this cell. Cells in the progeny are the daughters, granddaughters, etc., of the progenitor (or precursor) cell, which originate from division events (either self-renewal, asymmetric or symmetric) in any compartment over time, not including the progenitor cell itself. It is a summary statistic of the process which quantifies the proliferative potential of a single cell in compartment \(C_i\), and its offspring. For example, in Fig. 3, we represent a particular realisation of the stochastic process with \(G_1=8\).The mean number of cells in the progeny of a progenitor cell, \(m_i={\mathbb {E}}[G_i]\), for any initial compartment of interest \(i\in \{1,\ldots ,N\}\), can be obtained with first-step arguments by conditioning on the next event that occurs in the stochastic process. This approach leads to the following system of equations$$\begin{aligned} &\Delta _1m_1= {} \Lambda _1m_2 + 2(\lambda _1+a_1+s_1), \\& \Delta _im_i= {} \Lambda _i m_{i+1} +\xi _{i} m_{i-1} +2(\lambda _i+a_i+s_i),\quad i\in \{2,\ldots ,N-1\}, \\& \Delta _Nm_N= {} \xi _N m_{N-1} + 2\lambda _N. \end{aligned}$$The system above can be expressed in matrix form with the column vectors \(\textbf{m}=\left( m_1, \ldots , m_N \right) ^T\) and \(\textbf{b}=\left( 2(\lambda _1+a_1+s_1), \ldots ,\right.\) \(\left. 2(\lambda _{N-1}+a_{N-1}+s_{N-1}),\right.\) \(\left. 2\lambda _N\right) ^T\), as follows$$\begin{aligned} \textbf{J}\, \textbf{m} = \textbf{b}, \end{aligned}$$
(10)
with a tri-diagonal coefficient matrix$$\begin{aligned} \textbf{J}= \left( \begin{array}{cccccc} \Delta _1&{}- \Lambda _1 &{} 0 &{} 0 &{} \cdots &{} 0\\ – \xi _{2} &{} \Delta _2 &{} -\Lambda _2&{} 0&{} \cdots &{}0\\ 0 &{} -\xi _3 &{} \Delta _3&{} -\Lambda _3&{} \cdots &{}0\\ \vdots &{}\ddots &{}\ddots &{} \ddots &{}\ddots &{} \vdots \\ 0&{} \cdots &{} 0 &{}- \xi _{N-1} &{} \Delta _{N-1} &{} -\Lambda _{N-1}\\ 0 &{} \cdots &{} 0 &{}0 &{} – \xi _{N} &{} \Delta _{N} \end{array} \right) . \end{aligned}$$
(11)
Figure 3A realisation of the stochastic process tracking the progeny of a single progenitor cell which starts in compartment \(C_1\). The cell tracked (see Section “Lifeline analysis”) is shown as striped. Whenever the tracked cell divides, we follow one of the daughter cells selected at random (see division in \(C_3\)). For each cell, the colour indicates the compartment where it is at any given time. Here, the tracked cell dies in \(C_3\) (brown), while its progeny continues up to \(C_4\). In this example, \(G_1 = 8 = G_1(1)+G_1(2)+G_1(3)+G_1(4)=4+0+2+2\).One can exploit the tri-diagonal structure of \(\textbf{J}\) to obtain an explicit, or recursive, solution. In particular, by following a Gaussian forward-elimination backward-substitution approach, such as the Thomas algorithm47,48, one can obtain the recursive equations$$\begin{aligned} m_N= & {} \rho _N,\quad m_i \ =\ \rho _i – \gamma _i m_{i+1},\quad i\in \{1,\ldots ,N-1\}, \end{aligned}$$
(12)
where \(\gamma _1 = -\Delta ^{-1}_1 \Lambda _1\), \(\rho _1 = 2\Delta ^{-1}_1(\lambda _1+s_1+a_1)\), and$$\begin{aligned} \gamma _i= & {} -\frac{\Lambda _i}{\Delta _i+\xi _{i} \gamma _{i-1} },\quad i \in \{2,\ldots ,N-1\}, \\ \rho _i= & {} \frac{2(\lambda _i+s_i+a_i) + \xi _{i} \rho _{i-1}}{\Delta _i+\xi _{i} \gamma _{i-1}},\quad i\in \{2,\ldots ,N\}. \end{aligned}$$Using backward-substitution, this recursive scheme leads to the explicit solution$$\begin{aligned} m_i= & {} \sum _{j=i}^N\, (-1)^{j-i}\rho _j \, \left( \prod _{l=i}^{j-1} \gamma _l \right) ,\quad i\in \{1,\ldots ,N\}, \end{aligned}$$
(13)
where \(\prod _{l=i}^{i-1} \gamma _l = 1\). A condition on the parameters arises during the implementation of the recursive scheme$$\begin{aligned} \Delta _1>0,\quad \Delta _i+\xi _i\gamma _{i-1}> & {} 0,\quad i\in \{2,\ldots ,N\}, \end{aligned}$$so that the mean values \(m_1,\ldots ,m_N\) are finite and non-negative, for all \(i\in \{1,\ldots ,N\}\). This ensures that the number of cells in the progeny of a progenitor cell is finite with probability one, \({\mathbb {P}}(G_i<+\infty )=1\), since \(m_i={\mathbb {E}}[G_i]={\mathbb {E}}[G_i \vert G_i<+\infty ]{\mathbb {P}}(G_i<+\infty )+{\mathbb {E}}[G_i \vert G_i=+\infty ]{\mathbb {P}}(G_i=+\infty )\). In the irreversible case, where \(\xi _i=0\) for \(i\in \{2,\ldots ,N\}\), the solution above simplifies to$$\begin{aligned} m_i=2\sum ^{N}_{j=i} \frac{\lambda _j+a_j+s_j}{\Delta _{j}}\left( \prod _{l=i}^{j-1} \frac{\Lambda _l}{\Delta _l}\right) ,\quad i\in \{1,\ldots ,N\}, \end{aligned}$$
(14)
where for \(j=N\), we set \(a_N=s_N=0\). For \(i=N\), the empty product above is equal to one, so that \(m_N=\frac{2\lambda _N}{\mu _N-\lambda _N}\). For the irreversible model the condition to have finite and non-negative solutions becomes \(\Delta _i>0\) for all \(i\in \{1,\ldots ,N\}\). Direct inspection of Fig. 2 shows that the condition \(\Delta _i>0\) avoids unlimited accumulation of cells in compartment \(C_i\).Probability generating functionLet us now go beyond the mean number of cells. For the irreversible model we can consider the probability generating function of \(G_i\),$$\begin{aligned} \Phi _i(z)={\mathbb {E}}(z^{G_i})=\sum _{k=0}^{+\infty }{\mathbb {P}}(G_i=k) z^k. \end{aligned}$$The variable \(G_i\) counts the cells in the progeny of a progenitor cell starting in compartment \(C_i\), arising from division events (self-renewal, asymmetric and symmetric division), not including the progenitor cell itself. To include the progenitor cell, one can define \(S_i \equiv G_i+1\), and denote the new generating function by \(\Psi _i(z)={\mathbb {E}}(z^{S_i})\). The total expectation law over all possible first events implies that$$\begin{aligned} {\mathbb {E}}(z^{S_i}) = \sum _{E_j} {\mathbb {E}}(z^{S_i} \> \vert \> \text {event} \> E_j){\mathbb {P}}(E_j), \end{aligned}$$
(15)
with \(E_j\in \{\)death, differentiation, self-renewal, asymmetric division, symmetric division\(\}\). This leads to$$\begin{aligned} \Psi _i(z)= & {} \frac{\mu _i}{\Sigma _i}z + \frac{\nu _i}{\Sigma _i}z \Psi _{i+1}(z) + \frac{\lambda _i}{\Sigma _i} z (\Psi _i(z))^2 + \frac{a_i}{\Sigma _i}z \Psi _i(z) \Psi _{i+1}(z)+\frac{s_i}{\Sigma _i}z (\Psi _{i+1}(z))^2, \end{aligned}$$
(16)
with \(\Sigma _i=\mu _i+\nu _i+\lambda _i+a_i+s_i\). Since \(S_i=G_i+1\), one can write \(z \Phi _i(z) = \Psi _i(z),\) so that$$\begin{aligned} \lambda _i z^2 \Phi _{i}^2(z) + (a_i z^2 \Phi _{i+1}(z) -\Sigma _i) \Phi _{i}(z) + \Phi _{i+1}(z)(\nu _iz +s_i z^2 \Phi _{i+1}(z))+ \mu _i\ =\ 0. \end{aligned}$$
(17)
The probability generating functions above agree with the mean values obtained earlier. In particular, by differentiating with respect to z, and setting \(z=1\), we have$$\begin{aligned} (\Sigma _i-a_i-2\lambda _i) \Phi _{i}'(1) = 2(a_i+\lambda _i+s_i) + (\nu _i +a_i+2s_i) \Phi _{i+1}'(1), \end{aligned}$$
(18)
which can be solved recursively, leading to$$\begin{aligned} {\mathbb {E}}(G_i) = \Phi _{i}'(1)&= 2 \sum _{j=i}^N \frac{a_j+\lambda _j+s_j}{\Delta _j} \left( \prod _{l=i}^{j-1} \frac{\Lambda _l}{\Delta _l} \right) , \end{aligned}$$
(19)
in agreement with Eq. (14). We also note that when \(i=N\), the cells in the progeny of a single progenitor cell in compartment \(C_N\) arise from a linear birth-and-death process, which in discrete time has death probability, \(p_d=\frac{\mu _N}{\mu _N+\lambda _N}\), and birth probability, \(p_b=\frac{\lambda _N}{\mu _N+\lambda _N}\). Then, a first-step argument for the random variable \(S_N\) leads to$$\begin{aligned} \Psi _N(z)= & {} p_d z + p_b z \Psi _N^2(z) \ =\ z \phi _N(\Psi _N(z)), \end{aligned}$$with \(\phi _N(z)=p_d+ p_b z^2\). We can then write$$\begin{aligned} \Phi _N(z)=\frac{\Psi _N(z)}{z}= \frac{z \> \phi _N(\Psi _N(z))}{z} = \phi _N(z \> \Phi _N(z)), \end{aligned}$$which has solution$$\begin{aligned} \Phi _N(z) = p_d\frac{1-\sqrt{1-4x}}{2x}\quad \text {where}\quad x=p_dp_bz^2. \end{aligned}$$
(20)
We find \({\mathbb {P}}(G_N=k)\) as the coefficient of the power of \(z^k\), using the property of the Catalan numbers49. In particular, since$$\begin{aligned} {\sum _{n=0}^{\infty }C_nx^n = \frac{1-\sqrt{1-4x}}{2x}}, \end{aligned}$$we obtain$$\begin{aligned} \begin{aligned} {\Phi _N(z)}&{= p_d\sum _{n=0}^{\infty }C_n(p_dp_b)^nz^{2n}= \sum _{n=0}^{\infty }\frac{(2n)!}{(n+1)!n!}p_d^{n+1}p_b^nz^{2n}= p_d + z^2p_d^2p_b + z^42p_d^3p_b^2 + z^65p_d^3p_b^4 + \cdots .} \end{aligned} \end{aligned}$$
(21)
Note that \({\mathbb {P}}(G_N=k+2)/{\mathbb {P}}(G_N=k) = \frac{4k}{k+3}p_dp_b\). Assuming \(\lambda _N<\mu _N\), we verify that \({\mathbb {E}}(G_N) = \Phi _N'(1) = \frac{2\lambda _N}{\mu _N-\lambda _N}\).Compartmental analysis of the progenyCells in the progeny of a single progenitor can belong to different compartments, as shown in Fig. 3. For a progenitor cell starting in compartment \(C_i\), compartments \(C_j\), \(j\in \{1,2,\ldots ,N\}\) with greater proliferative potential will contribute more to \(G_i\). The proliferative potential of compartment j depends on the parameters \(\lambda _j,a_j,s_j,\mu _j,\nu _j,\xi _j\), but also on the number of cells arriving into that compartment. It is of interest then to write \(G_i=\sum _{j=1}^NG_i(j)\), with \(G_i(j)\) the number of cells in compartment \(C_j\) which belong to the progeny of the progenitor cell from compartment \(C_i\). For example, for the stochastic realisation of Fig. 3, \(G_1=G_1(1)+G_1(2)+G_1(3)+G_1(4)=4+0+2+2=8\).One can follow similar arguments to the ones in Section “The progeny of a single progenitor cell” to compute the mean quantities \(m_i(j)\equiv {\mathbb {E}}[G_i(j)]\). In particular, for an initial compartment \(C_i\), a first-step argument yields the following equations$$\begin{aligned} & \Delta _i m_{i}(i)= {} 2 \lambda _i + a_i +\Lambda _i\, m_{i+1}(i) +\xi _i m_{i-1}(i), \\& \Delta _im_{i}(i+1)= {} 2 s_i + a_i +\Lambda _i\, m_{i+1}(i+1) +\xi _{i}m_{i-1}(i+1), \\& \Delta _im_{i}(j)= {} \Lambda _im_{i+1}(j) +\xi _{i} m_{i-1}(j),\quad j\notin \{i,i+1\}, \end{aligned}$$where we implicitly set \(\xi _{1} = 0\), and \(\forall j \in \{1, 2, \ldots , N\}\), \(m_{N+1}(j)=m_N(N+1)=m_0(j) = 0\) for notational convenience. Making use of a recursive approach one can show that for any \(j\in \{1,\ldots ,N\}\), we have$$\begin{aligned} m_N(j)= & {} \rho _N(j),\quad m_i(j) \ =\ \rho _i(j) – \gamma _i m_{i+1}(j),\quad i\in \{1,\ldots ,N-1\}, \end{aligned}$$where \(\gamma _1 = -\Delta ^{-1}_1 \Lambda _1\), \(\rho _1(j) = \Delta ^{-1}_1 d_{(1,j)}\), and$$\begin{aligned} \gamma _i= & {} -\frac{\Lambda _i}{\Delta _i+\xi _{i} \gamma _{i-1} },\quad \quad i \in \{2,\ldots ,N-1\}, \\ \rho _i(j)= & {} \frac{d_{(i,j)}+ \xi _{i} \rho _{i-1}(j)}{\Delta _i+\xi _{i} \gamma _{i-1}},\quad i,j\in \{2,\ldots ,N\}, \end{aligned}$$with$$\begin{aligned} d_{(i,j)}= & {} {\left\{ \begin{array}{ll} 2\lambda _i + a_i, &{} \quad \text {if} \quad j=i, \\ 2s_i+a_i, &{} \quad \text {if} \quad j=i+1,\\ 0, &{} \quad \text {otherwise.} \end{array}\right. } \end{aligned}$$This recursive scheme leads to the solution$$\begin{aligned} m_{i}(j) =\sum _{k=i}^N(-1)^{k-i}\rho _ k(j)\, \left( \prod _{p=i}^{k-1} \gamma _p \right) , \end{aligned}$$
(22)
where \(\prod _{p=i}^{i-1} \gamma _p = 1\). This expression simplifies further for the irreversible model, since \(\xi _i=0\) for all \(i\in \{2,\ldots ,N\}\). In this instance, \(m_i(j)=0\) whenever \(i>j\), and we can write$$\begin{aligned} m_{i}(i)= & {} \Delta _i^{-1} (2 \lambda _i + a_i ),\\ m_{i}(j)= & {} \Delta ^{-1}_{j-1} \left( \prod _{p=i}^{j-2}\Delta ^{-1}_p \Lambda _p\right) \left( d_{(j-1,j)}+ d_{(j,j)}\Delta ^{-1}_{j} \Lambda _{j-1}\right) ,\quad j\in \{i+1, \ldots ,N\}, \end{aligned}$$for any \(i\in \{1,\ldots ,N-1\}\). Finally, we note that \(\prod _{p=i}^{i-1}\Delta ^{-1}_p \Lambda _p=1\), and \(m_{N}(N) = \frac{2\lambda _N}{\Delta _N}\).Lifeline analysisIn the previous sections we have analysed several summary statistics of the progeny of a progenitor cell. This section is devoted to summary statistics related to the history of a tracked cell, extending the analysis proposed in Ref.42, in the situation when all cell divisions are self-renewals; that is, \(s_i=a_i=0\) for all \(i\in \{1,\ldots ,N\}\). We construct lifelines containing a single cell at any time, starting in compartment \(C_i\). The line continues when that cell divides (by selecting one daughter cell at random) and terminates when the cell dies. In any given compartment, this cell can divide, move to another compartment, or die. Since the only type of division is self-renewal, we use the theoretical “trick” of identifying one of the daughter cells as a continuation of the mother, and the other as a new cell. The dynamics can then be represented by the CTMC \({{{\mathcal {Y}}}}=\{Y(t):\ t\ge 0\}\) over the state space \({{{\mathcal {S}}}}=\{ C_1, C_2,\ldots , C_N,\emptyset \}\), where Y(t) represents the state of the cell at time t; the cell is either in some compartment, \(C_j\), or it has died (state \(\emptyset\)). A schematic representation of the process \({{{\mathcal {Y}}}}\) is given in Fig. 4, and a particular realisation of this stochastic process can be seen in Fig. 3, where the tracked cell is shown with stripes.Figure 4Representation of the process \({{{\mathcal {Y}}}}\) to follow the fate of a single cell.We study the stochastic process with the following summary statistics:
The duration, \(T_i\), of the lifeline, starting in compartment \(C_i\), $$\begin{aligned} T_i= & {} \inf \{t\ge 0:\ Y(t)=\emptyset \vert Y(0)=C_i\}, \end{aligned}$$ which measures the survival potential of cells in the system depending on their initial compartment,
The number of divisions along a path through the family tree, \(D_i\), which quantifies the proliferative capacity of cells according to their initial compartment, and
The probability that the lifeline ends in a compartment \(C_j\); that is, \(\beta _i(j)={\mathbb {P}}(Y(T_i-\Delta t)=C_j)\) for small enough \(\Delta t\).
Lifespan of a tracked cellLet \(\tau _i={\mathbb {E}}[T_i]\) be the mean lifespan of a cell starting in compartment \(C_i\). By conditioning on the events the cell can undergo in the initial compartment, we obtain the following recursive equations$$\begin{aligned} &(\mu _1+\nu _1) \tau _{1}= {} \nu _1 \tau _{2} +1,\\& (\mu _i+\nu _i+\xi _i) \tau _{i}= {} \nu _i \tau _{i+1} + \xi _i \tau _{i-1} + 1,\quad i\in \{1,\ldots ,N-1\},\\& (\mu _N+\xi _N) \tau _{N}= {} \xi _N \tau _{N-1} + 1. \end{aligned}$$The previous equations have the general solution$$\begin{aligned} \tau _i =\sum _{k=i}^N(-1)^{k-i}{{\bar{\rho }}}_ k\, \left( \prod _{p=i}^{k-1} {{\bar{\gamma }}}_p \right) , \end{aligned}$$where \(\prod _{p=i}^{i-1} {{\bar{\gamma }}}_p = 1\). We have made use of the notation \({{\bar{\gamma }}}_1 = -{{\bar{\Delta }}}^{-1}_1 \nu _1\), \({{\bar{\rho }}}_1 = {{\bar{\Delta }}}^{-1}_1\), with$$\begin{aligned} {{\bar{\gamma }}}_i= & {} -\frac{\nu _i}{{{\bar{\Delta }}}_i+\xi _{i} {{\bar{\gamma }}}_{i-1} },\quad i \in \{2,\ldots ,N-1\}, \\ {{\bar{\rho }}}_i= & {} \frac{1+ \xi _{i} {{\bar{\rho }}}_{i-1}}{{{\bar{\Delta }}}_i+\xi _{i} {{\bar{\gamma }}}_{i-1}},\quad i\in \{2,\ldots ,N\}, \end{aligned}$$where \({{\bar{\Delta }}}_i = \mu _i + \nu _i +\xi _i\). A simpler solution is obtained in the irreversible model$$\begin{aligned} \tau _i =\sum _{k=i}^N \frac{1}{\mu _{k}+\nu _{k}} \, \left( \prod _{p=i}^{k-1} \frac{\nu _{p}}{\mu _{p}+\nu _p} \right) ,\quad i\in \{1,\ldots ,N\}, \end{aligned}$$where \(\prod _{p=i}^{i-1} \frac{\nu _{p}}{\mu _{p}+\nu _p} = 1\), and we set \(\nu _N=0\) for notational convenience.A similar approach allows one to compute the Laplace-Stieltjes transform of \(T_i\), and any of its higher order moments50. In particular, one could obtain a system of equations, via first-step arguments, for the Laplace-Stieltjes transform of \(T_i\), \(\phi _i(z)={\mathbb {E}}[e^{-zT_i}]\). The different order moments can be computed via successive differentiation of the transform, while the probability density function of \(T_i\) can be computed via its numerical inversion51. For example, in the irreversible model, the second-order moment of the lifespan of a cell starting in compartment \(C_i\) is given by$$\begin{aligned} {\mathbb {E}} [T_{i}^2] = \sum ^{N}_{j=i} R_j \left( \prod _{r=i}^{j-1}\frac{\nu _r}{\mu _r+\nu _r} \right) ,\quad i\in \{1,\ldots ,N-1\}, \end{aligned}$$where \(\prod _{p=i}^{i-1} \frac{\nu _r}{\mu _r+\nu _r} = 1\), \(R_i= \frac{2(\nu _i \tau _{i+1}+1)}{(\mu _i+\nu _i)^2 }\), and \(R_N = \frac{2}{\mu _N^2}\). We note that if the cell starts in the last compartment \(C_N\), \({\mathbb {E}} [T_N^2]=2\mu _N^{-2}\), since \(T_N\sim \text {Exp}(\mu _N)\).Number of divisions along a lifelineLet us define \(D_i\), the number of division events along the history of the tracked cell, starting in compartment \(C_i\). We compute the average value, \(\eta _i={\mathbb {E}}[D_i]\), with the first-step formula$$\begin{aligned} & (\lambda _1+\mu _1+\nu _1)\eta _1= {} \lambda _1(\eta _1+1)+\nu _1\eta _{2},\\& (\lambda _i+\mu _i+\nu _i+\xi _i)\eta _i= {} \lambda _i(\eta _i+1)+\nu _i\eta _{i+1}+\xi _i\eta _{i-1},\quad i\in \{1,\ldots ,N-1\},\\& (\lambda _N+\mu _N+\xi _N)\eta _N= {} \lambda _N(\eta _N+1)+\xi _N\eta _{N-1}. \end{aligned}$$The previous equations have solution$$\begin{aligned} \eta _i =\sum _{k=i}^N(-1)^{k-i}{{{\tilde{\rho }}}}_ k\, \left( \prod _{p=i}^{k-1} {{\bar{\gamma }}}_p \right) , \end{aligned}$$
(23)
where \(\prod _{p=i}^{i-1} {{\bar{\gamma }}}_p = 1\), \({{\tilde{\rho }}}_1 = {{\bar{\Delta }}}^{-1}_1 \lambda _1\), and$$\begin{aligned} {{\tilde{\rho }}}_i= & {} \frac{\lambda _i+ \xi _{i} {{\tilde{\rho }}}_{i-1}}{{{\bar{\Delta }}}_i+\xi _{i} {{\bar{\gamma }}}_{i-1}},\quad i\in \{2,\ldots ,N\}, \end{aligned}$$with \({{\bar{\Delta }}}_i\) and \({{\bar{\gamma }}}_i\) defined as above. In the irreversible model, this expression simplifies to$$\begin{aligned} {\eta _i }=\sum _{k=i}^N \frac{\lambda _k}{\mu _{k}+\nu _{k}} \, \left( \prod _{p=i}^{k-1} \frac{\nu _{p}}{\mu _{p}+\nu _p} \right) ,\quad i\in \{1,\ldots ,N-1\}, \end{aligned}$$where \(\prod _{p=i}^{i-1} \frac{\nu _{p}}{\mu _{p}+\nu _p} = 1\). We note that in this case \(\eta _N=\lambda _N\mu _N^{-1}\), since \(D_N\sim \text {Geometric} \left( \frac{\lambda _N}{\mu _N+\lambda _N}\right)\).A division event may occur at any instant along the lifeline of the cell, which may visit different compartments over time. Thus, one can determine the proliferative potential of the cell during its (eventual) visit to each compartment by considering \(D_i=\sum _{j=1}^ND_i(j)\), where \(D_i(j)\) is the number of divisions which occur exactly in compartment \(C_j\). Average values, \(\eta _i(j)\equiv {\mathbb {E}}[D_i(j)]\), can be computed (again) with first-step arguments. For instance, in the irreversible model, one has$$\begin{aligned} \eta _i(i)= & {} \frac{\lambda _i}{\mu _i + \nu _i}, \\ \eta _{i}(j)= & {} \frac{\lambda _j}{\mu _j+\nu _j}\left( \prod _{k=i}^{j-1} \frac{\nu _{k}}{\mu _k+\nu _k}\right) ,\quad j \ge i+1, \end{aligned}$$for \(i\in \{1, \ldots , N\}\) and \(\eta _N(N+1)=0\). We note that the expression above is consistent with the interpretation that, in the irreversible model, \(D_i(j)\sim \text {Geometric}\left( \frac{\lambda _j}{\lambda _j+\nu _j+\mu _j}\right)\) restricted to the arrival of the cell to compartment \(C_j\). In general, one can write$$\begin{aligned} \eta _i(j)= & {} {\mathbb {E}}[D_i(j)]={\mathbb {E}}[D_i(j)\ \vert \ {cell\;ever\;visits }~ C_j]\>{\mathbb {P}}(cell\;ever\;visits ~ C_j)+{\mathbb {E}}[D_i(j)\ \vert \ cell\;never\;visits ~ C_j]\>{\mathbb {P}}(cell\;never\;visits ~ C_j), \end{aligned}$$and, since \({\mathbb {E}}[D_i(j)\ \vert \ {cell\;never\;visits} ~ C_j]=0\), the quantity \(\eta _i(j)\) accounts for the probability of the cell not visiting this compartment.More generally, we consider the complete probability distribution of \(D_i\), defined by \(\omega _i(n) \equiv {\mathbb {P}}(D_i =n)\), the probability that a single cell starting in compartment \(C_i\) divides exactly n times before it dies or leaves the system, for any non-negative integer n. One can show that$$\begin{aligned} \omega _{i}(n)= & {} \sum ^{N}_{k=i} (-1)^{k-i}{{\hat{\rho }}}_k(n) \left( \prod _{j=i}^{k-1} {{\hat{\gamma }}}_j\right) , \quad i\in \{1,\ldots ,N\}, \,\, {n =0,1,2,3, \ldots }, \end{aligned}$$where \(\prod _{j=i}^{i-1} {{\hat{\gamma }}}_j= 1\). We have introduced the notation \({{\hat{\rho }}}_1(n) = (\lambda _1 \omega _1(n-1)+ \mu _1 1_{n=0}){{\hat{\Delta }}}_1^{-1}\), \({{\hat{\gamma }}}_1 = -\nu _1 {{\hat{\Delta }}}_1^{-1}\), and$$\begin{aligned} {{\hat{\rho }}}_i(n) = \frac{\lambda _i \omega _i(n-1) +\mu _i 1_{n=0}+ \xi _i {{\hat{\rho }}}_{i-1}(n)}{{{\hat{\Delta }}}_i+\xi _i {{\hat{\gamma }}}_{i-1}},\quad {{\hat{\gamma }}}_i = \frac{-\nu _i}{{{\hat{\Delta }}}_i+\xi _i {{\hat{\gamma }}}_{i-1}}, \end{aligned}$$with \({{\hat{\Delta }}}_i = \mu _i + \lambda _i + \xi _i + \nu _i\). We note that \(\omega _i(-1) = 0\) and that \(1_{n=0}\) is the indicator function, equal to 1 if \(n=0\), and zero otherwise. Thus, the probabilities \(\omega _i(n)\) can be computed recursively for increasing values of n, since \(\omega _i(n)\) depends on \(\omega _i(n-1)\).Cell deathThe probability that the tracked cell dies in compartment j is denoted by$$\begin{aligned} \beta _{i}(j)= & {} {\mathbb {P}}(Y(T_i-\Delta t)=C_j), \end{aligned}$$for small enough \(\Delta t\), and for any \(i,j\in \{1,\ldots ,N\}\). Once again, a first-step argument leads to the following recursive relations$$\begin{aligned}& (\mu _1+\nu _1)\beta _{1}(j)= \nu _1 \beta _{2}(j)+\mu _1 1_{j=1}, \\& (\mu _i+\nu _i+\xi _i)\beta _{i}(j)= \nu _i \beta _{i+1}(j)+\xi _i \beta _{i-1}(j)+\mu _i 1_{i=j},\\ &(\mu _N+\xi _N)\beta _{N}(j)= \xi _N \beta _{N-1}(j)+\mu _N 1_{j=N}, \end{aligned}$$with solution$$\begin{aligned} \beta _{i}(j)= & {} \sum ^{N}_{k=i} (-1)^{k-i} \bar{\rho }_k(j) \left( \prod _{p=i}^{k-1} \bar{\gamma }_p\right) \quad i,j\in \{1,\ldots ,N\}, \end{aligned}$$where \(\bar{\rho }_1(j) = \mu _1 \bar{\Delta }_1^{-1} 1_{j=1}\) for \(j\in \{1,\ldots ,N\}\), and$$\begin{aligned} \bar{\rho }_i(j) = \frac{\mu _i 1_{i=j} +\xi _i \bar{\rho }_{i-1}(j)}{\bar{\Delta }_i +\xi _i\bar{\gamma }_{i-1}},\quad i\in \{2, \ldots , N\},\ j\in \{1, \ldots ,N\}. \end{aligned}$$For the irreversible model, for \(i\in \{1,\ldots ,N\}\), we have$$\begin{aligned} \beta _{i}(i)= & {} \frac{\mu _i}{\mu _{i} + \nu _{i}},\quad \beta _{i}(j) \ =\ \frac{\mu _j}{\mu _j+\nu _j}\prod _{k=i}^{j-1} \frac{\nu _{k}}{\mu _k+\nu _k}, \quad j\in \{i+1, \ldots ,N\}. \end{aligned}$$We conclude this analysis with a comment. A particular advantage of obtaining analytical expressions for the summary statistics is that they allow one to explicitly compute sensitivities, \(\partial \beta _i(j)/\partial \theta\), or elasticities, \((\partial \beta _i(j)/\partial \theta )/(\beta _i(j)/\theta\)), with respect to model parameters of interest, \(\theta\). This can be rather useful when considering a complex model (with many parameters), as illustrated in Section ’Tracking a thymocyte during its development”. A local sensitivity analysis of this kind provides a quantification of the impact that small perturbations of model parameters can have on a given summary statistics of interest, and is particularly relevant when a subset of the model parameters are being estimated from experimental data sets, and thus, will carry inherent uncertainties. Finally, we refer the reader to the Supplementary Material for extra details on different calculations related to this section.