Previous studies (Cowan and Staudte, 1986; Sandler et al., 2015) found that simple inheritance rules, where interdivision times are correlated from one generation to another through a single parameter, cannot explain the lineage correlation patterns seen in experimental single-cell data. To address this issue, we propose a unified framework where the interdivision time is determined by a number of cell cycle factors that represent hidden variables such as cell cycle phase lengths, protein levels, cell growth rate or other unknowns (Figure 1c), that each have their own inheritance pattern.

The states of the cell cycle factors is assumed to be a vector ${\mathit{y}}_p={(y_{p,1},y_{p,1},\mathrm{\dots },y_{p,N})}^{\top }$ that determine the interdivision time of a cell with index $p$ via

Inheritance from mother to daughter of the $N$ cell cycle factors is described by a nonlinear stochastic Markov model on a lineage tree:

where $m\text{in}\mathbb{N}$ denotes the mother cell index and $2m$ and $2m+1$ the daughter cell indices. $f:{ℝ}_{+}^N\to {ℝ}_{+}$ and $\mathit{g}:{ℝ}_{+}^N\to {ℝ}_{+}^N$ are possibly nonlinear functions that model the dependence of the interdivision time on cell cycle factors and the inheritance process. ${\mathit{e}}_p={(e_{p,1},e_{p,2},\mathrm{\dots },e_{p,N})}^{\top }$ is a noise vector for which the pair ${\mathit{e}}_{2m},{\mathit{e}}_{2m+1}$ are identically distributed random vectors with covariance matrix independent of $m$. A non-zero covariance between these noise vectors can account for correlated noise of sister cells. We implicitly assume symmetric cell division such that the deterministic part of the inheritance dynamics $\mathit{g}$ is identical between the daughter cells. Note that we choose Equation 1a to be deterministic since division noise can be modelled by adding one more cell cycle factor that does not affect inheritance dynamics $\mathit{g}$.

The general model (Equation 1a and b) includes many known cell cycle models as a special case. For example, the interactions between cell cycle factors could model cell size control mechanisms (Appendix 1 - Section A6.1), the coordination of cell cycle phases (Appendix 1 - Section A6.3), or deterministic cues, such as periodic forcing of the cell cycle (Appendix 1 - Section A7.1), or coupling of the circadian clock to cell size control (Appendix 1 - Section A7.2).

The full model can only be solved for specific choices of $f$ and $\mathit{g}$, and these functions are generally unknown in inference problems. To overcome this limitation, we assume small fluctuations resulting in an approximate linear stochastic system (see Appendix 1 - Section A1 for a derivation) involving the interdivision time

The vector of cell cycle factor fluctuations ${\mathit{x}}_p={(x_{p,1},x_{p,2},\mathrm{\cdots },x_{p,N})}^{\top }$ obeys

Here, $\overline{\tau }$ is the stationary mean interdivision time, $\mathit{\theta }$ is the $N\times N$inheritance matrix and ${\mathit{z}}_{2m}$ and ${\mathit{z}}_{2m+1}$ are two noise vectors of length $N$ that capture the stochasticity of inheritance dynamics and differentiate the sister cells (Figure 2a). We denote the $N\times N$ covariance matrices ${\mathit{S}}_1=\mathrm{Var}({\mathit{z}}_{2m})=\mathrm{Var}({\mathit{z}}_{2m+1})$ and ${\mathit{S}}_2=\mathrm{Cov}({\mathit{z}}_{2m},{\mathit{z}}_{2m+1}),\text{for all}m\text{in}\mathbb{N}$ of the noise terms $\mathit{z}$ (and $\mathit{e}$) in individual cells and between sister cells, respectively. The noise terms are independent for all other family relations. The cell cycle factor fluctuations are scaled such that $\mathit{\alpha }={({\alpha }_1,{\alpha }_2,\mathrm{\dots },{\alpha }_N)}^{\top }$ is a binary vector of length $N$ made up of 1 s and 0 s depending on whether the function $f$ determining the interdivision time has dependence on a given cell cycle factor (see Appendix 1 - Section A1 for details). Under this scaling each cell cycle factor has a positive effect on the interdivision time, and hence we do not distinguish between factors with positive or negative effects on interdivision time.

Figure 2

Download asset Open asset

When the special case of a single cell cycle factor ($N=1$) is considered, the inheritance matrix model system reduces to a well-known model with correlated division times (Cowan and Staudte, 1986; Staudte et al., 1984; Staudte, 1992; Staudte et al., 1997), and we will refer to this case as simple inheritance rules (see also Appendix 1 - Section A5). In the following, we will explore the correlation patterns generated by multiple cell cycle factors.

Here, we define a correlation pattern to be the correlation coefficients of pairs of cells on a lineage tree. Here we introduce a function $\rho (k,l)$ which we call the generalised tree correlation function:

where ${\tau }_k$ and ${\tau }_l$ are the interdivision times of cells in the pair $(k,l)$, and $s_{\tau }$ is the interdivision time variance. The coordinate $(k,l)$ describes the distance in generations from each cell in the pair to their shared nearest common ancestor (Figure 2b and c). We have derived a closed-form formula for $\rho (k,l)$ (Equation M3) in Materials and Methods - ‘Analytical solution of the inheritance matrix model’; (see Appendix 1 - Section A3 for a full derivation) as a weighted sum of powers of the inheritance matrix eigenvalues $\lambda$:

with

We observe that the eigenvalues determine the dependence of the tree correlation function on $k$ and $l$, while the noise matrices ${\mathit{S}}_1$ and ${\mathit{S}}_2$ determine their relative weights $w_{ij}$ (see Equation 5).

Our theoretical analysis reveals three distinct correlation patterns that can be generated by the inheritance matrix model (further details in Materials and methods - ‘Analysis of tree correlation patterns’). These can be classified by the eigenvalues of the inheritance matrix $\mathit{\theta }$: (i) if the inheritance matrix exhibits real positive eigenvalues, we observe an aperiodic pattern (Figure 2d); (ii) if the inheritance matrix has real eigenvalues with at least one negative eigenvalue, we observe an alternator pattern (Figure 2e); and (iii) if there is a pair of complex eigenvalues we observe an oscillator pattern (Figure 2f). An intuitive interpretation of the eigenvalue decomposition is that it transforms the cell cycle factors into effective factors inherited independently. Hence, the inheritance matrix is diagonal in this basis. However, the analogy is limited to the case where the inheritance matrix is symmetric and the eigenvalues are real. For simplicity, we will focus on models with two cell cycle factors and note that in higher dimensions ($N\ge 3$), the correlation patterns involve a mixture of the three patterns discussed in detail in this section (Appendix 1—figure 6c, d and g,h).

To demonstrate the aperiodic correlation pattern, we utilise an inheritance matrix with positive real eigenvalues (Figure 2d). Characteristically, the modelled interdivision time correlations decay to zero as the distance to the most recent ancestor increases (Figure 2g) since the eigenvalues in Equation 4 are bounded between 0 and 1. To look more closely at the patterns on the tree, we utilise two reductions of the generalised tree correlation function. These are the lineage correlation function ($\rho (k,l)$ for $k$ or $l=0$) and the cross-branch correlation function ($\rho (k,l)$ for $k=l$). We look at these functions for continuous $k,l$ to visualise better the patterns that occur down the lineage and across the branches of the tree. The lineage correlation function gives the correlation dynamics as you go down the lineage tree, whereas the cross-branch correlation function gives the correlation dynamics as you move across neighbouring branches of the lineage tree. We observe that the interdivision time correlations decrease as we move both across generations and branches (Figure 2j).

In contrast, the alternator pattern generates oscillations with a fixed period of two generations in the lineage correlation function. The behaviour is typically observed for cell cycle factors with negative mother-daughter correlations (Appendix 1 - Section A6.1). In this case, we have at least one negative eigenvalue and thus Equation 4 will alternate between positive and negative values for successive generations, producing the period two oscillation. We demonstrate this correlation pattern for the generalised tree correlation function (Figure 2h) using a diagonal $\mathit{\theta }$ matrix (Figure 2e). We observe alternating correlations across generations in the lineage correlation function, and the continuous interpolation of the cross-branch correlation function (Figure 2k). Although the period is fixed to two generations, the amplitude of the correlation oscillation varies with the absolute magnitude of the eigenvalues (Materials and methods - ‘Analysis of tree correlation patterns’).

To investigate the oscillator correlation pattern, we propose a hypothetical inheritance matrix $\mathit{\theta }$ with eigenvalues $\mathit{\lambda }=(D{e}^{+i\frac{2\pi }{P}},D{e}^{-i\frac{2\pi }{P}})$ which are complex for $D,P\ne 0$ and $P\ne \frac{2}{k},k\text{in}\mathbb{Z}$ (Figure 2f). The parameters $P$ and $D$ control the period and the respective damping of an underlying oscillator, i.e., the limit $D\to 1$ leads to an undamped oscillation and $D\to 0$ corresponds to an overdamped oscillation (see Materials and methods - ‘Determining the period of correlation oscillations from the eigenvalues’ for details). Correspondingly, the graph of the generalised tree correlation function (Figure 2i) shows clear oscillations across generations. These correlation oscillations are also evident in the lineage correlation function but are absent in the cross-branch correlation function (Figure 2l). However, oscillations are possible in the cross branch correlation function for other choices of $\mathit{\theta }$ with complex eigenvalues (see model fits in Figure 3 and Methods - ‘Analysis of tree correlation patterns‘). In summary, the qualitative behaviour of the interdivision time correlation patterns can be studied using the eigenvalue decomposition of the inheritance matrix $\mathit{\theta }$.

Our analysis shows that of the three specified patterns, only the oscillator pattern cannot arise from simple inheritance rules. This is because it requires at least two inherited cell cycle factors ($N\ge 2$) for the inheritance matrix to possess complex eigenvalues. We therefore asked whether the oscillator pattern is necessary for the cousin-mother inequality to be satisfied. We find that this is not the case, but instead, all three correlation patterns can be compatible with the cousin-mother inequality if $N\ge 2$. To demonstrate this, we choose three specific two-dimensional inheritance matrices $\mathit{\theta }$ that produce the required eigenvalue structure (Figure 2d–f). We then use these matrices with our analytical solution for the generalised tree correlation function (Materials and methods - ‘Analytical solution of the inheritance matrix model’) to map the regions where the cousin-mother inequality can be satisfied (Figure 2m–o). Interestingly, we find that oscillations can arise even in parameter regions that violate the cousin-mother inequality (Figure 2o). We conclude that both the cousin-mother inequality and the oscillator pattern are sufficient but not necessary conditions to rule out simple inheritance rules.

To understand which datasets can be explained by simple inheritance rules, we fit the one-dimensional model ($N=1$) to six publicly available lineage tree datasets (Appendix 1—table 1) using Bayesian methods (Materials and methods - ‘Data analysis and Bayesian inference of the inheritance matrix model’). These datasets were chosen as they each had a sufficient number of cells for correlation analysis and covered a broad range of cell types. We found that the model fit is poor for the datasets that display the cousin-mother inequality, which is the case for cyanobacteria, clock-deleted cyanobacteria, neuroblastoma and human colorectal cancer cells (Appendix 1—figure 1a–f). Despite not obeying the cousin-mother inequality, the fit is also poor for mouse embryonic fibroblasts (Appendix 1—figure 1f) as the median inferred correlation lies outside the 95% confidence intervals for both the grandmother and cousin correlations which are included in the model fit, and the confidence intervals for the data vs the credible intervals from the inference show minimal overlap (Appendix 1—figure 2f). Another inequality may be violated in this dataset that cannot be explained using the one-dimensional model, suggesting that the absence of the cousin-mother inequality cannot rule out more complex division rules. The only cell type that has a good fit for the one-dimensional model is mycobacteria (Appendix 1—figure 1c). We thus conclude that the majority of the datasets must be described by higher dimensional inheritance dynamics of multiple cell cycle factors.

We asked whether the correlation patterns are better described by a two-dimensional inheritance matrix model. Bayesian inference (Materials and methods - ‘Data analysis and Bayesian inference of the inheritance matrix model’) produced a good model fit for all six datasets (Figure 3a–f) for the two-factor inheritance matrix model, within relatively narrow error bars of mother, grandmother, sister and cousin correlations (Appendix 1—table 1). The credible intervals from the Bayesian inference matched the confidence intervals of correlations used for fitting (Appendix 1—figure 2). We quantified the quality of our fits using the Akaike information criterion (AIC) (Materials and methods - ‘Data analysis and Bayesian inference of the inheritance matrix model’, (Equation M11)) for each dataset and compared these to the one-dimensional model (Appendix 1—table 1). The AIC estimates the goodness of fit with a penalty for model complexity allowing us to select the simplest model that explains the data. The AIC values indicate that the inheritance matrix model with two cell cycle factors provides the simplest fit for all cell types used here, except for the mycobacteria data where simple inheritance rules provided an equally good fit with a significant reduction in the number of model parameters. We expected the AIC to select the two dimensional model where the cousin-mother inequality was satisfied such as in cyanobacteria, clock-deleted cyanobacteria, neuroblastoma, and human colorectal cancer cells. The match with the two-factor inheritance matrix model in fibroblasts was less obvious.

Figure 3

Download asset Open asset

Crucially, we find that the model has a good predictive capacity for correlations further down the lineage tree. For each pattern, we show several samples from the conditional posterior distribution (solid and shaded lines) to illustrate fits of the lineage correlation and cross-branch correlation function (Figure 3a–f). For all datasets except neuroblastoma, the curves also intercept the great-grandmother and great-great-grandmother correlations that were not used for fitting (Figure 3a–d and f), and bootstrapped confidence intervals from the data overlapped with the credible intervals obtained from Bayesian inference (Appendix 1—figure 2). We then asked which correlation patterns underlie the data. To assess this, we calculated the eigenvalues of each posterior sample of the inheritance matrix to categorise the aperiodic, alternator and oscillator patterns (Figure 3a–f, bar charts). We found that in every dataset, the dominant correlation pattern was identifiable with probabilities well above 50%, except for mycobacteria (Figure 3c) that was better described by simple inheritance rules (Appendix 1—figure 1c).

Cyanobacteria, (Figure 3a), human colorectal cancer (Figure 3d), and mouse embryonic fibroblasts (Figure 3f) display a dominant oscillator pattern, but we see that their lineage correlation functions exhibit widely different periodicities. For example, the posterior lineage correlation for cyanobacteria displays a higher frequency oscillation than those in human colorectal cancer cells and fibroblasts. Clock-deleted cyanobacteria (Figure 3b) and mycobacteria (Figure 3c) display a dominant alternator pattern which could be induced by strong sister correlations. We see that clock-deleted cyanobacteria (Figure 3b) has a 100% alternator pattern in contrast to the 100% oscillator pattern seen for wild type cyanobacteria, suggesting that the deletion of the clock gene has completely transformed the correlation pattern and has abolished the underlying oscillation. Neuroblastoma (Figure 3e) displays a dominant aperiodic pattern. The predictive capacity for this cell type is weaker than for the other datasets, which we assume is due to the tight confidence interval in the correlations. Despite this discrepancy, we find that the inheritance matrix model produces excellent fits and has good predictive capacity for all other cell types studied in this work.

We next ask which mechanisms are responsible for generating the observed correlation patterns. The Bayesian inference used for model fitting (Materials and methods - ‘Data analysis and Bayesian inference of the inheritance matrix model’) samples parameters using a MCMC Gibbs sampler. The Gibbs sampler can be thought of as a random walk in parameter space that settles around parameter regions with high likelihood. We found that the explorations of the Gibbs sampler did not settle in a particular parameter subspace but meandered off to explore vast areas of the parameter space without improving the likelihood values (Appendix 1—figure 3a and b). Such behaviour is expected when model parameters are not identifiable and the posterior distribution of parameters cannot be efficiently sampled (Rannala, 2002; Raue et al., 2013).

To provide further evidence of unidentifiablity, we obtained four histograms of a single parameter of the inheritance matrix for different initialisations. The four distributions are very different (Figure 4a), showing that the random walk does not settle to a stationary distribution. We further observe that the mean squared displacement increases without bound (Figure 4b) showing that the sampling does not settle in a particular subset of the parameter space. In contrast to the individual parameters, the sampled posterior distribution of the eigenvalues is consistent across the averages (Figure 4c) and their mean squared displacement converges rapidly (Figure 4d). We note that unidentifiability arises for the inheritance matrix model with multiple cell cycle factors and does not feature for simple inheritance rules (Appendix 1 - Section A5). This ultimately demonstrates that the interdivision time correlation patterns do not identify a single set of inheritance parameters, but rather need to be described by a distribution of inheritance mechanisms.

Figure 4

Download asset Open asset

Clock-deleted cyanobacteria and neuroblastoma both satisfy the cousin-mother inequality (Figure 1b), which indicates that at least two cell cycle factors are responsible for the corresponding correlation patterns. The eigenvalues of the inheritance matrix concentrate in different regions of the admissible parameter space (Figure 4e), suggesting the correlation patterns that generate the cousin-mother inequality are distinct. For the clock-deleted cyanobacteria dataset, we found that all posterior samples were consistent with an alternator correlation pattern, while most posterior samples presented aperiodic correlation patterns in neuroblastoma (Figure 3b and e bar charts).

We hypothesised that different inheritance models generate these patterns. To verify this hypothesis and since we cannot identify the cell cycle factors directly, we computed the mother-daughter correlations between the two hidden cell cycle factors. Since the order of factors is interchangeable, we only distinguish between mother-daughter correlations between the same ($\text{corr}(x_{m,i},x_{2m,i})$ and $\text{corr}(x_{m,i},x_{2m+1,i})$ for $i=1,2$) and alternate factors ($\text{corr}(x_{m,i},x_{2m+1,j})$ and $\text{corr}(x_{m,i},x_{2m,j})$ for $i\ne j=1,2$). The resulting posterior distributions revealed distinct correlation patterns of cell cycle factor correlations for clock-deleted cyanobacteria and neuroblastoma (Figure 4f). For clock-deleted cyanobacteria, we predict that at least one factor has a negative mother-daughter correlation while its cross-correlation with the other factor must be positive; while the correlations are of opposite sign for neuroblastoma (Figure 4f). We sketch influence diagrams that summarise these relationships between factors (Figure 4g and h). Thus, the different interdivision time correlation patterns observed for clock-deleted cyanobacteria and neuroblastoma stem from distinct hidden correlation patterns of cell cycle factor fluctuations.

We observe that the lineage correlation functions of cyanobacteria, human colorectal cancer cells, and fibroblasts exhibit vastly different correlation oscillation periods (Figure 3). Next, we are interested to see whether the oscillations seen in these datasets are compatible with biological oscillators known to affect cell cycle control.

Correlation oscillations and underlying rhythms can exhibit vastly different periods

The period of the correlation oscillation is related to the location of the eigenvalues of the inheritance matrix on the complex plane. We consider an eigenvalue $\lambda$ of the inheritance matrix. In terms of the mean interdivision time $\overline{\tau }$, the correlation period ${\displaystyle T_0}$ is:

(6) ${\displaystyle T_0=\overline{\tau }\frac{2\pi }{|\text{Arg}(\lambda )|}\ge 2\overline{\tau },}$

and the inequality means that the period ${\displaystyle T_0}$ is always greater than twice the mean interdivision time $\overline{\tau }$. More generally, there is an oscillation period associated with each eigenmode of the inheritance matrix, but the period is infinite for real eigenvalues, and thus only complex eigenvalues generate correlation oscillations. This inequality follows from Equation 6 using $|\text{Arg}(\lambda )|\le \pi$. However, known biological oscillators that influence cell cycle control often have periods less than twice the mean interdivision time, such as stress response regulators (Harper et al., 2018; Stewart-Ornstein et al., 2017) and gene expression oscillations (William et al., 2007; Gao et al., 2014; Whitfield et al., 2002). How can relatively slow observed correlation oscillations be compatible with much faster biological oscillators underlying the cell cycle?

The resolution to this issue is that the period of the correlation oscillation does not always match the frequency of the underlying oscillator. Instead there are a number of possible oscillator periods $T_n$ compatible with the correlation oscillation period ${\displaystyle T_0}$ (Appendix 1 - Section A4) given by:

(7) ${\displaystyle T_n=\frac{\overline{\tau }{T}_0}{|\overline{\tau }+n{T}_0|},}$

for $n\text{in}\mathbb{Z}$. This phenomenon, that the same correlation oscillation can be explained by multiple underlying oscillators, can be understood using the intuition in Figure 5a.

Figure 5

Download asset Open asset

Code available at https://github.com/pthomaslab/Lineage-tree-correlation-pattern-inference (copy archived at swh:1:rev:dc69bbce5ce909813d7d4356c9fd2da045e02c79; Hughes, 2022).

From Equation 2b and $\mathrm{IE}[{\mathit{z}}_p]=0,\text{for all}p\text{in}\mathbb{N}$, we see that the vector of cell cycle factors has zero mean $\mathrm{IE}[{\mathit{x}}_p]=0$. Its $N\times N$ covariance matrix $\mathbf{\sum }=\mathrm{Cov}({\mathit{x}}_p,{\mathit{x}}_p)$ satisfies a discrete-time Lyapunov equation:

From the solution of Equation M1, we compute the variance of the interdivision time

and the generalised tree correlation function $\rho (k,l)$ (see Appendix 1 - Section A3 for a detailed derivation) given by:

where $\mathit{\omega }(k,l)={\mathit{\theta }}^k\mathbf{\sum }{\left({\mathit{\theta }}^l\right)}^{\top }+{\delta }_{k\ge 1}{\delta }_{l\ge 1}{\mathit{\theta }}^{k-1}{\mathit{S}}_2{\left({\mathit{\theta }}^{l-1}\right)}^{\top }$ with

To ensure that the lineage tree correlation pattern is stationary, we require ${\displaystyle \text{SR}(\mathit{\theta })<1}$ where $\text{SR}(\mathit{\theta })=\text{max}({\lambda }_1,{\lambda }_2,\mathrm{\dots },{\lambda }_N)$ is the spectral radius of $\mathit{\theta }$. This also ensures that the solutions to Equation M1; $\mathbf{\sum }$, ${\mathit{S}}_1$ and the function Equation M3 are unique and independent of the initial conditions.

The patterns of the generalised tree correlation function can be characterised through its eigendecomposition. The general decomposition proceeds through finding the matrix of eigenvectors $U$ of $\mathit{\theta }$ such that

is the diagonal matrix of eigenvalues. Defining ${\widehat{\mathit{S}}}_{1,2}=U{\mathit{S}}_{1,2}{U}^{\top }$ and $\widehat{\mathit{\alpha }}={(U^{-1})}^{\top }\mathit{\alpha }$, the solution to Equation M1 is given by

This result can then be used to find an explicit expression for the generalised tree correlation function:

where

Equation M7 can be rewritten as a superposition of patterns Equation 4 with weights given by Equation 5.

The pattern of the tree correlation function is thus governed by the eigenvalues of the inheritance matrix $\mathit{\theta }$: (i) if one eigenvalue, say ${\lambda }_1$, is positive then the factor ${\widehat{\omega }}_{11}(k,0)={\widehat{\omega }}_{11}(0,k)\propto {\lambda }_1^k$ contributing to lineage correlation decays monotonically. The factor ${\widehat{\omega }}_{11}(k,k)\propto {\lambda }_1^{2k}$ contributing to the cross-branch correlation decays twice as fast; (ii) if there is a negative eigenvalue, the factor ${\widehat{\omega }}_{11}(k,0)={\widehat{\omega }}_{11}(0,k)\propto {(-1)}^k{|{\lambda }_1|}^k$ alternates between negative and positive values with an envelope of ${|{\lambda }_1|}^k$, while the corresponding contribution to the cross-branch correlation decays monotonically with rate as ${|{\lambda }_1|}^{2k}$. Finally, if we have a pair of complex eigenvalues ${\lambda }_1={\lambda }_2^{\ast }=D{e}^{i\mathrm{\Omega }}$ then the factors ${\widehat{\omega }}_{i,j}(k,0)={\widehat{\omega }}_{i,j}(0,k)$ contributing to the lineage correlation function display damped oscillations with frequency $\mathrm{\Omega }$ and envelope $D^k$, while the factor ${\widehat{\omega }}_{12}(k,k)={\widehat{\omega }}_{12}^{\ast }(k,k)\propto D^{2k}$ and the factor ${\widehat{\omega }}_{11}(k,k)={\widehat{\omega }}_{22}^{\ast }(k,k)\propto D^{2k}{e}^{i2\mathrm{\Omega }k}$ oscillate with frequency $2\mathrm{\Omega }$.

We consider the case where the inheritance matrix $\mathit{\theta }$ has a pair of complex conjugate eigenvalues ${\lambda }^{\pm }=D{e}^{\pm i2\pi /P}$. The lineage correlation function then oscillates whenever $D\ne \{0,1\}$ and $P\ne \frac{2}{k},k\text{in}\mathbb{Z}$. The period of correlation oscillations per generation is given by

where ${\displaystyle \text{Arg}(\lambda )\text{in}(-\pi ,\pi ]}$ is the argument of the eigenvalue and $\ln (\cdot )$ is the complex logarithm. The former is the angle made between the line joining the origin and the eigenvalue $\lambda$ on the complex plane with the real axis. This means that $T_0/\overline{\tau }=P$ if and only if ${\displaystyle P>2}$. Otherwise, ${\displaystyle T_0}$ is calculated in terms of $P$ by Equation M9, (Appendix 1—figure 8).

We determined all pairs of cells in a lineage tree, sorted them by family relations $(k,l)$ and calculated the sample correlation coefficient of interdivision times (Equation 3). To maximise the number of samples used to calculate these correlations, an individual cell can appear in more than one pair. For example, if a cell had two cousins, it would be counted in two separate cousin pairs in the cousin-cousin correlation coefficient calculation. For training, we focus on the sample statistics $\widehat{\mathit{X}}=({\widehat{s}}_{\tau },{\{{\widehat{\rho }}_{(k,l)}\}}_{(k,l)\text{in}C})$ with $C=\{(1,0),(2,0),(1,1),(2,2)\}$ comprised of the interdivision time sample variance and four interdivision time sample correlation coefficients given by the mother-daughter, grandmother-granddaughter, sister-sister and cousin-cousin relations (Figure 2a). Note that ${\widehat{s}}_{\tau }$ is computed across all interdivision times used to calculate the correlation coefficients in each dataset. Errors are estimated using bootstrapping by re-sampling cell pairs with replacement 10,000 times. The resulting variances and correlation coefficients are given in Appendix 1—table 1.

The vector of inferred model parameters for the two-dimensional model is $\mathbf{\Theta }=(\mathit{\theta },{\mathit{S}}_1)$, where we fix $\mathit{\alpha }={(1,1)}^{\top }$ and ${\mathit{S}}_2=0$ for simplicity. A different choice of $\mathit{\alpha }$ did not affect our results (Appendix 1—figure 4). Since ${\mathit{S}}_1$ is symmetric, it consists of the $N$ variances and $N(N-1)/2$ correlation coefficients between the components of $\mathit{z}$. Thus for $N=2$ the inheritance matrix model has seven free parameters to be estimated. We assumed that the log-likelihood for these statistics is the sum of square errors:

which is equivalent to assuming that the sample variance and correlation coefficients are normally distributed for large sample sizes. We calculate the interdivision time variance $s_{\tau }$ and the generalised tree correlation function $\rho (k,l)$ from Equation M2 and Equation M3. Note that Equation M2 is the interdivision time variance from a tree where all lineages have the same number of generations, which approximates the variance across all cells in the observed trees (Appendix 1—table 3). For simplicity, we neglected possible correlations between the sample statistics in $\widehat{\mathit{X}}$ and used bootstrapped estimates for the standard deviation of the sample statistics ${\widehat{\sigma }}_{{\widehat{s}}_{\tau }}$ and ${\widehat{\sigma }}_{{\widehat{\rho }}_{(k,l)}}$ (Appendix 1—table 1). Note that the likelihood is independent of the mean since it is irrelevant for the correlation pattern. We assumed a flat prior with support restricted to ${\displaystyle \text{SR}(\mathit{\theta })<1}$ and ${\mathit{S}}_1$ positive semi-definite to guarantee the existence of a stationary correlation pattern.

The numerical implementation uses the adaptive Gibbs‐sampler implemented in the Julia library Mamba.jl (Smith, 2018). For each dataset, we sample 11 million parameter sets which include a burn-in transient of 1 million samples. These samples are removed before analysis of the output.

For model comparison we use the AIC (Akaike, 1974) given by

where $k$ is the number of model parameters and $\ln \left(\widehat{ℒ}\right)$ is the maximum value of the log-likelihood function given by Equation M10. For ${\mathit{S}}_2\ne 0$, the inheritance matrix model has $k=d(1+2d)$ parameters where $d$ is the number of cell cycle factors in the model. For ${\mathit{S}}_2=0$ the number of parameters reduces to $k=\frac{1}{2}d(1+3d)$.

Here, we will derive the inheritance matrix model given by Equation 2a, and b in the main text. We assume that the fluctuations in the hidden cell factor dynamics are small, which leads to a computationally efficient approximation.

Firstly, in the limit of zero fluctuations, all cells must be identical. Hence, all cell cycle factors are equal to their means ${\displaystyle \mathit{\mu }=({\mu }_1,{\mu }_2,\cdots ,{\mu }_N{)}^{\mathrm{\top }}=\mathbf{E}(y_p)}$ and similarly for the noise vectors ${\displaystyle \beta =({\beta }_1,{\beta }_2,\cdots ,{\beta }_N{)}^{\mathrm{\top }}=\mathbf{E}(e_{2m})=\mathbf{E}(e_{2m+1})}$ in Equation 1b. From Equation 1a, and b we then find that

which can be efficiently solved for ${\displaystyle \overline{\tau }}$ and ${\displaystyle \mu }$ using standard numerical methods.

Secondly, we can decompose the interdivision time and the cell cycle factor vector into their respective mean and fluctuating components by

Denoting the index of the present cell by p and the one of its mother by m, we can expand f and g around the limit of zero fluctuations and we obtain to leading order

where

Using this expansion and Equation A2 in Equation 1a and b of the main text we arrive at

where we have set ${\mathit{e}}_p=\mathit{\beta }+{\tilde{\mathit{z}}}_p$ and ${\tilde{\mathit{z}}}_p={({\tilde{z}}_{p,1},{\tilde{z}}_{p,2},\mathrm{\cdots },{\tilde{z}}_{p,N})}^{\top }$ giving the fluctuations around the mean for the noise vectors. Comparing Equation A7 with Equation A1 and collecting terms to leading order, we obtain the linearised system:

Next, we define the diagonal scaling matrix $\mathbf{\Gamma }$ with non-zero elements as

for $i=1,2,\mathrm{\dots },N$. Using the rescaled noise sources $\mathit{z}=\mathbf{\Gamma }\tilde{\mathit{z}}$, we find the rescaled inheritance matrix $\mathit{\theta }$ and $\mathit{\alpha }$-coefficients

The rescaled cell cycle factor fluctuations ${\mathit{x}}_m=\mathbf{\Gamma }{\tilde{\mathit{x}}}_m$ follow Equation 2 of the main text and we reach rescaled variance-covariance matrices ${\mathit{S}}_1$ and ${\mathit{S}}_2$ as follows

Here we analyse the effect of nonlinearity on the interdivision time correlation patterns. For simplicity we consider a single cell cycle factor and follow the same lines as in Appendix 1 - Section A1, Equation A9, while including terms of order $x^2$. This leads to the expansion interdivision time and factor fluctuations

where $\tilde{\theta }$ is the Jacobian of the cell cycle factor dynamics, as before, and $\tilde{H}=g^{\prime \prime }({\mu }_x)$ is the Hessian. From the second equation we obtain

Defining ${\tilde{\mathit{X}}}_p=({\tilde{x}}_p,{\tilde{x}}_p^2)$ and ${\tilde{\mathit{Z}}}_p=({\tilde{z}}_p,{\tilde{z}}_p^2)$, combining Equation A14 and Equation A16 and rescaling variables as in Equation A12 and Equation A13, we find the extended inheritance matrix model

where

Here $H=\tilde{H}\tilde{\alpha }/\tilde{\beta }$ and $\beta =1$ if $\tilde{\beta }\ne 0$, and analogously, $H=\tilde{H}\tilde{\alpha }$ and $\beta =0$ if $\tilde{\beta }=0$. Hence, the interdivision time correlation patterns with small to moderate fluctuations can be described through an extended linear system (Equation A7) that includes nonlinear terms $x_p^2$. These additional terms can be interpreted as cell cycle factors forming binary complexes. The presence of these complexes increases the number of cell cycle factors and extends the eigenvalue spectrum of the effective inheritance matrix $\mathrm{\Theta }$ by ${\theta }^2$. Hence, the presence of complexes leads to mixed correlation patterns. For example, for a single cell cycle factor, the eigenvalues of $\mathbf{\Theta }$ are $(\theta ,{\theta }^2)$, which corresponds to an alternator pattern for ${\displaystyle \theta <0}$. More generally, we may expect that nonlinear patterns can be described through mixtures of aperiodic, alternator, and oscillatory patterns. For example, the complex eigenvalue spectrum of an oscillator pattern ($e^{\pm i2\pi /P}$) will include powers of complex eigenvalues ($e^{\pm i4\pi /P}$) resulting in harmonics of the fundamental correlation oscillation frequency similar to higher order harmonics observed in single-cell time-series of the circadian clock (Martins et al., 2016; Thomas et al., 2013).

In this section we derive an analytical expression for the generalised tree correlation function. This gives the Pearson correlation coefficient in interdivision time for any pair of related cells. We start with the equation for the Pearson correlation coefficient, and from there derive a formula for the interdivision time covariance using the known properties of the cell cycle factors $\mathit{x}$. From this, we can derive the general formula for the correlation coefficient between any related cell pair.

We associate a cell pair with an index $(k,l)$ which measures the distance to the nearest common ancestor as given in ‘The inheritance matrix model reveals threedistinct interdivision time correlation patterns’ (Figure 2a). From this, we denote their interdivision time fluctuations as ${\tau }_k^{\prime }$ and ${\tau }_l^{\prime }$ respectively. The Pearson correlation coefficient between these fluctuations is given by

where $s_{{\tau }^{\prime }}$ is the variance of the interdivision time fluctuations.

The interdivision time fluctuations ${\tau }_k^{\prime }$ and ${\tau }_l^{\prime }$ are calculated from the vector of rescaled cell cycle factor fluctuations ${\mathit{x}}_k$ as given in ‘A general inheritance matrix model provides a unified framework for lineage tree correlation patterns’, giving the equations

Substituting Equation A20 into Equation A19, we obtain a formula for in terms of the cell cycle factor fluctuations and the coefficients alone

Since ${\mathit{x}}_k$ and ${\mathit{x}}_l$ are identically distributed in steady state, we have that $\mathrm{Var}({\mathit{x}}_k)=\mathrm{Var}({\mathit{x}}_l)=\mathrm{Cov}(\mathit{x},\mathit{x})=\mathbf{\sum }$ as specified in Materials and methods - ‘Analytical solution of the inheritance matrix model’‘. We can write $\rho (k,l)$ now as

where ${\mathit{\alpha }}^{\top }\mathbf{\sum }\mathit{\alpha }$ gives the variance of the interdivision time fluctuations ${\tau }^{\prime }$.

Using the model Equation 2 we can write the formula for the $\mathit{x}$ vectors for the two cells in the cell pair $(k,l)$ as

where cells $k$ and $l$ have mother cells $k-1$ and $l-1$ respectively. The two cells are sisters if and only if their subscripts are both equal to 1, meaning they share a mother cell. Using recurrence of the model, we can write these equations as

where ${\mathit{x}}_0$ is the vector of cell cycle factors for the most recent common ancestor for a cell pair given by $(k,l)$.

All that remains is to derive a function for $\mathrm{Cov}({\mathit{x}}_k,{\mathit{x}}_l)$ which we will denote $\mathit{\omega }(k,l)$. We calculate $\mathit{\omega }(k,l)$ as follows using expectations:

where ${\mathit{\mu }}_{{\mathit{x}}_k}$ and ${\mathit{\mu }}_{{\mathit{x}}_l}$ are the mean vectors of ${\mathit{x}}_k$ and ${\mathit{x}}_l$ respectively which are both equal to 0, giving

To find $\mathrm{IE}\left[{\mathit{x}}_k{\mathit{x}}_l^{\top }\right]$ in terms of the model parameters, we substitute in Equation A26 for ${\mathit{x}}_k$ and ${\mathit{x}}_l$ and get

The noise term fluctuations $\mathit{z}$ are only correlated if the cells are sisters, which only occurs when the distance we have $(k,l)=(1,1)$. So for the summations above, we exclude all terms except where $i=j=1$. Doing this and expanding we get

where ${\delta }_{k\ge 1}$ and ${\delta }_{l\ge 1}$ are given in Equation M4. We also have that

The matrix $\mathrm{Cov}({\mathit{x}}_0,{\mathit{x}}_0)$ is equivalent to the covariance matrix for any $\mathit{x}$, giving $\mathrm{Cov}({\mathit{x}}_0,{\mathit{x}}_0)=\mathbf{\sum }$. This gives

Similarly we have,

As $\mathrm{IE}(\mathit{z})=0$, and $\mathrm{Cov}({\mathit{z}}_{2m},{\mathit{z}}_{2m+1})={\mathit{S}}_2$ as stated in Materials and methods - ‘Analytical solution of the inheritance matrix model’, we obtain,

Equation A31 therefore becomes:

Substituting Equation A37 back into Equation A29 we get

giving us the final equation for $\mathit{\omega }(k,l)$. Using the above equation in Equation A23, we obtain Equation M3 of the Methods.

The period of correlation oscillation as observed in the lineage correlation functions is given by Thomas et al., 2018. We can reveal the underlying oscillator periods by shifting the inferred period ${\displaystyle T_0}$ to obtain a smaller period $T_n$. This means that shorter periods would produce the same inferred period in the lineage correlation function when sampled at the original frequency of once per cell cycle (Figure 5a).

The oscillator periods are obtained by adding or subtracting multiples of $2\pi$ to the argument of the eigenvalue which results in the new argument being in the same position in the complex plane. The oscillator period $T_n$ with shift $n\text{in}\mathbb{Z}$ is therefore given by

Taking Equation A40 and substituting in Equation 6, we obtain $T_n$ in terms of ${\displaystyle T_0}$ as Equation 7.

We consider the limiting case of a single cell cycle factor ($N=1$) resulting in simple inheritance rules. This situation could model a growth factor that can either increase or decrease interdivision times of cells depending on the monotonicity of $f$ in Equation 1. The analytical solution (Equation M7) of the inheritance matrix model then reduces to

where $w=1+{\delta }_{k\ge 1}{\delta }_{l\ge 1}\frac{S_2}{S_1}\frac{(1-{\theta }^2)}{{\theta }^2}$. First, we observe that, given single cell measurements of the mother-daughter correlation coefficient $\rho (0,1)$, the daughter-daughter correlation coefficient $\rho (1,1)$ and the variance $s_{\tau }$, the parameters $\theta$, S₁ and S₂ are uniquely identifiable:

Thus measurements of the variance, lineage- and cross-branch correlations fully determine the parameters. The tree correlation function is, however, independent of $f$, which means that the interdivision time correlation pattern carries no information whether the growth factor increases or decreases growth. The reason for this indifference is that cell cycle factors are identified only by their fluctuation pattern, i.e., for each cell cycle factor whose fluctuations increase interdivision time $x$, we could define another cell cycle factor fluctuation that decrease interdivision time $-x$. We accounted for this unidentifiability issue trough a similarity transformation using the scaling matrix $\mathbf{\Gamma }$ in Equation A12 and Equation A13 that transforms all cell cycle factor fluctuations to increase interdivision time. Of course, this unidentifiabiliy could be removed through explicitly measuring the involved cell cycle factors.

To further investigate the output of the inheritance matrix model, we propose multiple models of known cell cycle control mechanisms, and map them to our inheritance matrix model framework. All cell size models assume symmetric division.