Introduction

At each division, an initial cell often referred to as the mother, first duplicates its genetic material and cellular components, then it halves into two daughter cells. Depending on how the internal components are partitioned, we refer to symmetric division when the daughters are identical to the mother, while we speak of asymmetric cell division if the two daughter cells differ in size, cellular components, and/or fate. Emerging experimental observation suggests that asymmetrical partitioning plays an important role in cell-to-cell variability, cell fate determination, cellular aging, and rejuvenation1,2,3. For instance, Katajisto et al.4 have shown that mammalian epithelial stem-like immortalized cells partition mitochondria asymmetrically. Tripathi et al.5 use a computational model to show that partitioning is one of the main noise sources behind cancer cell plasticity and employs it to model epithelial-mesenchymal transition (EMT). In their model, partitioning noise increases the heterogeneity of the tumor environment and changes its ability to resist treatment5. In this respect, several theoretical and experimental works have shown that in populations of bacteria or cancer cells facing environmental changes, variability increases the probability that some individuals may survive the stress produced by a sudden change of the environment, e.g., the one produced by antibiotics or cancer treatments6,7,8,9,10,11.

Moreover, asymmetric cell division and differential centrosome inheritance have been found in human neuroblastoma cells12, while Hwang et al. showed that polarity factors control the segregation of subcellular vesicles during the division of colorectal cancer stem cells13. Similarly, mitochondria and endosomes are known to exhibit asymmetrical partitioning in yeast14,15,16,17,18,19 and in asymmetrically dividing cells, like mammalian epithelial stem-like immortalized cells20. This asymmetrical mitochondria partitioning enables cells to protect themselves from aging because it protects the cell progeny from accumulating misfolded proteins engulfed in the mitochondria4,18,19,21. Indeed, Bergmiller et al.22 found that biased binomial segregation helps bacterial populations to face antibiotic treatments, suggesting partitioning as a drug resistance enhancer. In addition, biased segregation of membrane proteins in bacteria has been linked to the formation of protein rafts which migrate differently to the older pole in the mother cell with respect to the novel one in the daughter cell, thus influencing the aging of the overall population23.

Whether similar mechanisms are also present in homogeneous non-differentiating cancer cell populations is still under scrutiny. In fact, Semrau et al.24 and Dunn et al.25 show how, in response to the state of their gene activation network, pluripotent cells differentiate in progenitor and terminal cells by gradually changing the abundance of their RNA and protein compounds. Several modeling approaches explained such processes assuming that cells are complex systems in a quasi-stable state25,26,27,28,29. In this framework, cell differentiation is the effect of sufficiently large perturbations of cell compounds that bring the cell from a quasi-stable state into a new one associated with a different cell type30. This perspective may be considered as a mathematical complex system formulation of the Waddington epigenetic landscape31. Nevertheless, given this general framework, it is unclear at the moment how a homogeneous non-differentiating cancer cell line can maintain its quasi-stable state under the effects of perturbations caused by partitioning asymmetry. It is also unclear if/why a homogeneous cell line would preserve a certain degree of asymmetrical partitioning.

Despite its biological importance, the quantification of asymmetry in cellular components at division is often based on microscopy time-lapse approaches that are time-consuming and not applicable in some systems. Furthermore, the division can be tracked only for few generations32 and often for only one compound at a time. To overcome these limitations, we developed a procedure to follow the proliferation of Jurkat cells, which is based on the measurement of how different cellular compounds partition. To do so, we combine a multicolor flow cytometry33,34 protocol with a computational procedure to identify cells belonging to different generations. Flow cytometry is a non-disruptive technique that enables the simultaneous acquisition of thousands of cells per second. The main idea of a multicolor flow cytometry approach is to have an initial population with two or more cell compartments stained with a live cell marker that does not affect the cell cycle and the partitioning process. Being a fast and non-invasive technique, flow cytometry has been successfully applied to study the growth and proliferation of very different cell types35,36.

Our method takes advantage of the high-throughput capability of flow cytometry to improve measurement accuracy. In particular, we analyzed the partition of cytoplasm, cellular membrane lipids, and mitochondria with the use of succinimidyl dye (CTV, CellTrace Violet), lipophilic dye (PKH, PKH Linker), and mitochondria dye (MITO, MitoTracker), respectively. To explain the flow cytometry data, we derived a statistical method based on the asymmetrical binomial partitioning model, which takes into account the possibility of having biases (and thus asymmetry) in the organelle partitioning process. In particular, comparing the fluorescence mean intensities and their fluctuations against the predictions of our model, we found that mitochondria and membrane lipids partition asymmetrically, while cytoplasm partitions symmetrically.

Thanks to our multi-label method, we are also able to investigate the correlation between the organelle partitioning. In fact, cellular partitioning involves several molecular compounds with different sizes and spatial distributions. In this perspective, the unexplained non-trivial correlations experimentally observed at the phenotypic level37,38,39,40 may be understood by the correlation between partitioning organelles and cell sub-components.

Results

Measuring organelle partitioning

To quantitatively measure the segregation of the cellular sub-components, we devised a multicolor flow cytometry experiment. The opportunity to stain more than one compartment in live cells allows us to follow how these compartments are partitioned in time and to study the guiding dynamical processes. In our experimental setup, we used three fluorescent dyes: Cell Tracer, PKH, and MitoTracker, which respectively label the cytoplasm, the cell membrane, and the mitochondria. See Fig. 1a for a schematic representation of the experimental setup, and “Methods” for more details. We note that Métivier et al.41 used MitoTracker Green FM to assess mitochondria content in cells. Indeed, this dye was described in flow cytometry protocols to establish mitochondria amount in human lymphoblastic cell lines42. Here, instead, we use MitoTracker in a sorting experiment to follow mitochondria partitioning in mammalian cell progeny.

Fig. 1: Experimental protocol.
figure 1

a Schematic representation of a fluorescence-activated cell sorter (FACS) working principle. Stained cells pass one at a time in front of a laser source that interacts with the fluorophores of the markers. The forward scattered (fsc), side scattered (ssc) light, together with the one emitted by the dye fluorophores is collected and analyzed. Eventually, a sorter divides cells according to certain thresholds on the measured intensities. b The dot plot shows how cells are collected on the basis of the highest expression level of PKH (PKH Linker, membrane lipid marker) and CTV (CellTrace Violet, cytoplasm marker). The ‘pre sorted’ panel presents a pre-sorting sample with the established collection gate and the ‘post sorted’ panel shows a post sorted sample. In this case, the sorted cells are 11.9% of the presorted population. c Sorted cells are then put in culture for subsequent acquisition of aliquots at the flow cytometer at different time points. d The histogram shows the CTV fluorescence intensity distribution that starting from the highest level at time zero (T0, red) decreases in time as each division occurs (orange, green, and blue). The gray-filled distribution represents the background intensity. e Two-marker dot plot. Each dot represents the measured intensity of PKH and CTV markers in a single cell. Points are colored from red to blue as the local density of points in the plot increases. Gray dashed lines mark the plane bisector. The top and left panels display the intensity marginal probability distributions of PKH and CTV intensities, respectively. f Analyzing the data with the expectation–maximization and Gaussian mixture algorithms enable us to identify different populations/generations (colored points). After initial guessing of the parameters, the algorithm optimizes the weight associated to each point. The position of the best-fit Gaussians is displayed by the concentric circles. g Same as in (e) but with cells clustered by generations according to the Gaussian mixture algorithm. Experimental data are represented as gray dots, while the position of the best-fit Gaussians is displayed by the concentric circles. Different colors correspond to different generations. Histograms represent the marginal distributions for the two markers.

Full size image

After staining the population with the fluorescence dyes, we selected an initial population with a narrowly peaked fluorescence distribution. This is an important step of the procedure. In fact, as all the cells of each generation divide into two daughter cells, the fluorescence distribution decreases its mean of one-half, and changes its variance as a consequence of the partitioning events; so the narrower the initial population the easier it is to discriminate the subsequent generations.

In fact, when the variance of the initial distribution is too broad the peaks of the subsequent generations tend to merge into a single peak. This effect becomes more important as the asymmetricity of the partitioning increases. The plots in Supplementary Fig. 1 show the simulated fluorescence distribution for two binomial partitioning models with different asymmetry parameters, p = 0.5 (low) and p = 0.42 (high). For higher asymmetry, the different generation peaks merge together even at smaller initial variances. See Supplementary Note 1 and 2 for further details.

To generate a narrowly peaked fluorescence distribution, we used fluorescence-activated cell sorting (FACS) that allows isolation of a given cell population based on the fluorescence signal of a specific color channel, subtracted by the background signal. Once established a gate of interest, based on each signal’s fluorescence intensity upon acquisition, all the cells enclosed in that gate are collected in vital and functional conditions. Cell sorter isolation ensures highly purified and sterile cells that can still be grown in culture for up to three weeks.

In this respect, we applied our protocol on human Jurkat T cells, that are a cancer cell line commonly used for proliferation studies. Jurkat cells have been verified especially for multiple dye staining and labeling efficiency, showing a low toxic response to dyes43. Indeed, in each sorting experiment we performed, cells were stained for viability, and were found highly viable (see the first plot of Supplementary Fig. 4 for the live cells gating strategy by dead-stain exclusion). Importantly, cells were passaged only twice upon defrosting to avoid excessive exposure to culture conditions prior to the experiment.

Figure 1b displays how cells are sorted around the mode of the PKH and CTV fluorescence distribution. Then, the sorted cells are put in culture (Fig. 1c). In many proliferation analyses, there is a tendency to synchronize cells using serum starvation to avoid that cells are in different stages of the cell cycle or temporally separated by a minimum of one division round. On the other hand, it has been found that Jurkat cells can be resistant to serum deprivation, and consequently that synchronization could be a source of error44,45. For this reason, we decided to not use whole-culture synchronized cells since we wanted to analyze cells in their proliferation state without perturbing the partitioning process experimentally. Furthermore, a recent work by Cooper46 discusses the issue of cell population synchronization, arguing that whole-culture synchronization methods are unreliable because they stop the cell cycle at a certain point, and they inhibit cell division. This procedure leads to biases in cell phenotype, for example, whole-culture synchronized cells exhibit broader cell size distribution in comparison to the unperturbed cell at the same cell-cycle point. Narrow sorting of stained cells may instead be effectively used to obtain unbiased cell synchronization.

Selection of a narrow initial distribution

Narrow fluorescent sorting indirectly produces cell synchronization. This is confirmed by two fundamental observable characteristics: (i) cells exhibit a narrow-size distribution, and (ii) cell populations show clear peaks corresponding to successive generations emerging and disappearing in correspondence of the inter-division time46. Indeed, because we sort the cells with narrow gates, centered on the bulk of the distribution of the fluorescence signal for each dye, we obtain a pool of cells that carry similar morphological characteristics. But the issue is more complex. We must take into consideration that a proliferating cell population is a dynamical system, consequently, cells in this population go through a sequence of phases characterized by specific morphology. Thus when we sort a cell population we have one of two situations: either we sort a homogeneous population of cells and consequently, we tend to synchronize them in a specific phase, or we sort an initially heterogeneous cell population and we risk selecting a particular more or less synchronized subset of cells. In our case, Jurkat cells are known to be a very homogeneous cell line especially under the specific condition set in the protocol. To further reinforce these considerations in Supplementary Note 2, we discuss the advantages of synchronizing a homogeneous population with cell sorting. In each experiment, we sampled cells in different initial conditions of phases and morphologies. As we obtain the same values of asymmetric partitioning (as quantified by the p parameter), we can conclude that these do not bias our result.

In our experiment, from the forward and side scattering data, we observe a narrow-size distribution compared to the unsorted cells size distribution. Furthermore, even if we do not sample several time points we find that after each successive inter-division time period, a new peak (associated with a new generation) emerges and contains the near-total majority of the cells of that generation (Supplementary Fig. 2). These observations confirm that the populations were indirectly synchronized by the narrow gate on the fluorescence intensity. However, it is important to point out that our method does not require cell-cycle synchronization but only a narrow sorting of the marker fluorescence distributions: the synchronization is a result of the adopted sorting procedure, but it is not required to measure asymmetry with our method.

Isolation of the different generations via Gaussian mixture

In practice, once the initial population is characterized, we start the measurement of the population dynamics by sampling a small portion of the growing population at an established kinetic. Then we estimate the fluorescence distribution of the population from the sampled distribution (see “Methods”). This enables us to follow how the fluorescence decreases at each analyzed time point (Fig. 1d). The simultaneous measurement of the fluorescence intensity for all the dyes enables us to associate a 2- or 3-dimensional point for each sampled cell, f = (f, g, h), where each dimension is the intensity of one of the fluorescent markers used. The fluorescent sub-components are divided into the new daughter cells within each new generation of cells; thus, at each division the average fluorescence of the new population will be half of the previous. Consequently, if we plot the base-two logarithm of f on a scatter plot, and if the initial variance of the fluorescence intensity of the cells at the first generation fg=0 is made sufficiently narrow (see Fig. 1b), we find that the different generations of cells form clusters on the scatter plot, and peaks on the histograms of the fluorescence intensities. Moreover, if we plot the histogram of the base-two logarithm of either f, g, or h, we see that each peak is at a unit distance from the next, and that these peaks have an approximate Gaussian shape, i.e., it is lognormal distributed in the linear scale (see Supplementary Fig. 5). In conclusion, if the initial variance of fg=0 is sufficiently small, Var[f0] ≪ E[f0], we get a sequence of peaks one for each generation g as shown in Fig. 1d.

We note that our procedure does not allow us to directly count the number of sub-components for each cell (mg), but it allows us to measure the fluorescent intensity (fg) that is detected by the flow cytometer. Where fg is the integral of the fluorescence signal collected over the time that the cell is exposed to the laser beam while it flows through the flow cytometer microfluidic channel. Moreover, we are implicitly assuming that the mean and variance are respectively proportional to E[fg] ~ aE[mg] + b and Var[fg] ~ a2Var[mg]; where a is the amount of intensity per unit sub-component, and b is the background fluorescence intensity.

To observe the change of variance due to the number of particles on the measured data we would have to produce a very narrow initial distribution that must be sufficiently small to allow us to observe the sub-generation peaks formed by the asymmetric segregation. Unfortunately, we are not able to get such a narrow peak due to the fact that its width results both from intrinsic and extrinsic factors43, respectively, the instrumental limit of the flow cytometer’s gate width, and the biological/chemical aspects of the given fluorescence dye. At each measurement, we get a distribution of fluorescence intensities, p(f), that is a mixture of the distributions corresponding on different generations:

$$p({{{{{{{\bf{f}}}}}}}})=\mathop{\sum}\limits_{{{{{{{{\bf{g}}}}}}}}}{\pi }_{{{{{{{{\bf{g}}}}}}}}}{{{{{{{\rm{P}}}}}}}}({{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\bf{g}}}}}}}}}).$$
(1)

Each peak would have an expected value such \({{{{{{{\mathrm{log}}}}}}}\,}_{2}({{{{{{{\rm{E}}}}}}}}[{{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\bf{g}}}}}}}}+{{{{{{{\bf{1}}}}}}}}}])={{{{{{{\mathrm{log}}}}}}}\,}_{{{{{{{{\bf{2}}}}}}}}}({{{{{{{\rm{E}}}}}}}}[{{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\bf{g}}}}}}}}}])-1\). Assume that \({{{{{{{\rm{P}}}}}}}}({{{{{{{\mathrm{log}}}}}}}\,}_{2}{{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\bf{g}}}}}}}}})\) is approximately Gaussian (Supplementary Fig. 5 shows how peaks can be approximated by Gaussian distributions in log scale), we can model the data as a Gaussian mixture

$$p({{{{{{{\mathrm{log}}}}}}}\,}_{2}{{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\bf{i}}}}}}}}})=\mathop{\sum}\limits_{{{{{{{{\bf{g}}}}}}}}}{\pi }_{{{{{{{{\bf{g}}}}}}}}}{{{{{{{\rm{N}}}}}}}}({{{{{{{\mathrm{log}}}}}}}\,}_{{{{{{{{\bf{2}}}}}}}}}{{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\bf{g}}}}}}}}}| {\tilde{\mu }}_{{{{{{{{\bf{g}}}}}}}}},{\tilde{\sigma }}_{{{{{{{{\bf{g}}}}}}}}})$$
(2)

where \({{{{{{{\rm{N}}}}}}}}({{{{{{{\mathrm{log}}}}}}}\,}_{2}{{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\bf{g}}}}}}}}}| {\tilde{\mu }}_{{{{{{{{\bf{g}}}}}}}}},{\tilde{\sigma }}_{{{{{{{{\bf{g}}}}}}}}})\) is a Gaussian distribution with mean \({\tilde{\mu }}_{g}\) and variance \({\tilde{\sigma }}_{g}\) (see “Methods”).

Figure 1e–g shows the Gaussian mixture protocol that assigns each point in the clouds to the corresponding generation. By performing measurements at regular intervals during the proliferation, we get sequences of samples that produce a snapshot of the population distribution at different times. Figure 2a shows an example of the population proliferation, followed through the intensities of the CTV and PKH markers. In this way, the sampled cells produce a series of clouds composed of data points, each of which represents a different generation. Each dot of these clouds corresponds to a cell that was part of the sample. The data points in Fig. 2a are presorted as routinely done in flow cytometry using the forward and side scatter to remove cell doublets and debris. One can spot the succession of peaks associated with the different generations that succeed one after the other in the population that grows and proliferates over time.

Fig. 2: Data analysis.
figure 2

a Dot plot of the intensity of PKH (PKH Linker, membrane lipid marker) vs CTV (CellTrace Violet, cytoplasm marker), measured at different times (hours) in one of the three experiments with simultaneous PKH and CTV markers. Points are colored from red to blue as the local density of points in the plot increases. Grey dashed lines mark the plane bisector. b Gaussian mixture fit from the estimate model parameters obtained from (a). Experimental data are represented as gray dots, while the position of the best-fit Gaussians are displayed by the concentric circles. Different colors correspond to different generations. The percentages, f, of cells in each generation, g, are reported in the insert as a bar plot. c Mean value of the intensity of the fluorescent markers as a function of the generations, obtained through the Gaussian mixture fitting. Blue triangles and red stars correspond to the mean of the intensity for cytoplasm and membrane dyes, respectively. The standard error of the mean is either represented with bars or is smaller than the size of the markers. d Same as in (c) but with the square root of the intensity variance. e Pearson correlations between the PKH and CTV intensities as a function of the generations. Green circles represent the mean correlation for each generation. The standard error of the mean is either represented with bars or is smaller than the size of the markers.

Full size image

Note that our model does not assume that cells divide all together and have the same cycle length, but instead it takes into account the possibility that cells have different cycle lengths. Nevertheless, it requires that each generation fluorescence peak must be detected by the Gaussian mixture algorithm, and that each peak has a measurable variance that is sufficiently constant in time. However, Supplementary Fig. 2 shows that different experiments presented consistent dynamics of each generation peak successions, indicating compatible cell cycles among populations.

Applying the Gaussian mixture fitting procedure with the expectation–maximization algorithm to each snapshot (see Fig. 2b), we can assign each point to one generation and then measure the mean (μg) and the variance (\({\sigma }_{g}^{2}\)) for each generation. As we can see from Fig. 2c, d, both the mean intensity and its fluctuations decrease as a function of the generations for the two displayed markers. The same behavior is shown with the mitochondria marker (see Supplementary Fig. 6). On the other hand, the Pearson correlation between the two intensities appears to remain low, but different from zero, during the proliferation, as shown in Fig. 2e.

Binomial partitioning model

To acquire quantitative insights on the partitioning mechanisms responsible for the observed dispersion of fluorescence intensities, we used an asymmetric binomial partitioning model. This model is a stochastic model that enables us to derive the parameters of the partitioning processes from experimentally measured means and variances of the intensity peaks. In the asymmetric binomial partitioning model, we assume that cell sub-components segregate randomly according to a Bernoulli process. To uniquely label all the cells in a progeny derived from a single progenitor cell with label 1, we consider that each cell i divides into two daughter cells that we can label unambiguously as 2i and 2i + 1. After g generations, the cell lineage generated from cell 1 forms a lineage tree. If we start with one cell, at generation g = 0, each generation is composed of ng = 2g cells (Fig. 3a).

Fig. 3: Theoretical predictions for the Bernoulli partition process.
figure 3

a Sketch of the division process of a cell having two kinds of sub-components marked with blue and red dyes. When a cell i with mi sub-components divides, m2i sub-components go to one daughter cell with probability p, and m2i+1 = mi − m2i sub-components go to the other daughter cell, with probability (1 − p). When the two daughters divide in turn, the sub-components inherited from the mother are further divided into niece cells. b Number of markers 1 and 2 per cell. Each cloud of points corresponds to a different generation. Points are colored from red to blue as the local density of points in the plot increases. Histograms represent the marginal distributions for the two markers. c Marginal distribution of marker 1 copy-numbers for generation 0 (dark red), 1 (red), and 2 (orange). For each generation, the dotted curves are associated with the different daughter cells subpopulations that originated from the previous generation, depending on the partition probability, p. Each sub-population can be associated with a product, \({{{{{{{\mathcal{F}}}}}}}}\), of fraction probabilities which depends on the division path. d Mean number of marked compounds per cell, μg divided by the mean at generation zero, μ, as a function of the generation, g for different possible values of the partition fraction, p. Solutions of Eq. (17) are represented with lines while dots stand for simulated values. e Square root variance of the number of marked compounds per cell, σg, divided by the square root variance at generation zero, μ, as a function of the generation, g for different possible values of the partition fraction, p.

Full size image

Given that the mother cell, i has mi sub-components that are stained with a fluorescent dye, the daughter cells 2i and 2i + 1 will have respectively m2i and m2i+1 (with mi = m2i + m2i+1) fluorescent sub-components. If we assume that segregation is symmetric and noiseless, the cell sub-components would always be divided exactly in half, i.e., m2i = m2i+1 = mi/2. Otherwise, if we assume that each component selects cell 2i or 2i + 1 with probability p, we find that m2i and m2i+1 are random variables and that, given mi, the distribution of m2i is binomial

$${{{{{{{\rm{P}}}}}}}}({m}_{2i}| {m}_{i},p)=\left(\begin{array}{c}{{m}_{i}}\\ {{m}_{2i}}\end{array}\right){p}^{{m}_{2i}}{(1-p)}^{{m}_{i}-{m}_{2i}}$$
(3)

with mean p ⋅ mi and variance p(1 − p)mi as depicted in Fig. 3a. In our flow cytometry experiment, at generation zero, we do not have a single mother cell but a population of cells characterized by a distribution over possible values of m1, P(m1), which has a specific average μ and variance σ2. The probability distributions of the inherited fluorescent sub-components in the first generation of daughter cells is given by the superposition of two subpopulations formed by the daughter cells 2 and 3 with, respectively, P(m2) and P(m3),

$${{{{{{{\rm{P}}}}}}}}({m}_{2})=\int {{{{{{{\rm{P}}}}}}}}({m}_{2}| {m}_{1},p){{{{{{{\rm{P}}}}}}}}({m}_{1})d{m}_{1}.$$

and

$${{{{{{{\rm{P}}}}}}}}({m}_{3})=\int {{{{{{{\rm{P}}}}}}}}({m}_{3}| {m}_{1},1-p){{{{{{{\rm{P}}}}}}}}({m}_{1})d{m}_{1}.$$

In general, given a cell i, and given the number of fluorescent sub-components P(mi), the distributions of the daughter cells 2i and 2i + 1 are

$${{{{{{{\rm{P}}}}}}}}({m}_{2i})=\int {{{{{{{\rm{P}}}}}}}}({m}_{2i}| {m}_{i},p){{{{{{{\rm{P}}}}}}}}({m}_{i})d{m}_{i},$$

and

$${{{{{{{\rm{P}}}}}}}}({m}_{2i+1})=\int {{{{{{{\rm{P}}}}}}}}({m}_{2i+1}| {m}_{i},1-p){{{{{{{\rm{P}}}}}}}}({m}_{i})d{m}_{i}.$$

Symmetric case

To derive a closed-form equation for the expectation and the variance of P(mi), we start with the simpler symmetric partition case (\(p=\frac{1}{2}\)), where P(m2i∣mi) and P(m2i+1∣mi) have average \({{{{{{{\rm{E}}}}}}}}[{m}_{2i}| {m}_{i}]={{{{{{{\rm{E}}}}}}}}[{m}_{2i+1}| {m}_{i}]=\frac{{m}_{i}}{2}\) and variance \({{{{{\mathrm{Var}}}}}}[{m}_{2i}| {m}_{i}]={{{{{\mathrm{Var}}}}}}[{m}_{2i+1}| {m}_{i}]=\frac{{m}_{i}}{4}\).

We begin with the expression of the expectation using the total expectation law, which is just a consequence of the double-integrals computed on the expectation value. Thus, setting μ1 = μ,

$${\mu }_{2i}={\mu }_{2i+1} ={{{{{{{\rm{E}}}}}}}}[{m}_{2i}]= \, {{{{{{{\rm{E}}}}}}}}[{{{{{{{\rm{E}}}}}}}}[{m}_{2i}| {m}_{i}]]=\int d{m}_{i}\int d{m}_{2i}{m}_{2i}{{{{{{{\rm{P}}}}}}}}({m}_{2i}| {m}_{i})\\ = \, E\left[\frac{{m}_{i}}{2}\right]=\frac{{\mu }_{i}}{2}$$
(4)

For the variance, the computation is a little more complex, but analogously, from the law of total variance, which is itself a consequence of total expectation law, and given σ1 = σ,

$${\sigma }_{2i}^{2} = \, {\sigma }_{2i+1}^{2}={{{{{{{\rm{Var}}}}}}}}[{m}_{2i}]={{{{{{{\rm{E}}}}}}}}[{{{{{{{\rm{Var}}}}}}}}[{m}_{2i}| {m}_{i}]]+{{{{{{{\rm{Var}}}}}}}}[{{{{{{{\rm{E}}}}}}}}[{m}_{2i}| {m}_{i}]]\\ = \, \frac{1}{4}{{{{{{{\rm{E}}}}}}}}[{m}_{i}]+\frac{1}{4}{{{{{{{\rm{Var}}}}}}}}[{m}_{i}]$$
(5)

For i = 1,

$${\sigma }_{2}^{2}={\sigma }_{3}^{2}={{{{{{{\rm{Var}}}}}}}}[{m}_{2}]=\frac{1}{4}{{{{{{{\rm{E}}}}}}}}[{m}_{1}]+\frac{1}{4}{{{{{{{\rm{Var}}}}}}}}[{m}_{1}]=\frac{\mu }{4}+\frac{{\sigma }^{2}}{4}$$
(6)

and for i = 1, 2

$${\sigma }_{4}^{2}={\sigma }_{5}^{2}={\sigma }_{6}^{2}={\sigma }_{7}^{2}=\frac{1}{4}{{{{{{{\rm{E}}}}}}}}[{m}_{2}]+\frac{1}{4}{{{{{{{\rm{Var}}}}}}}}[{m}_{2}]=\frac{3\mu }{16}+\frac{{\sigma }^{2}}{16}$$
(7)

Indeed, solving the iterative Eqs. (4) and (5), we get, for any cell i from generation g such that 2g ≥i > 2g+1,

$${\mu }_{i}={{{{{{{\rm{E}}}}}}}}[{m}_{i}]=\frac{\mu }{{2}^{g}}.$$
(8)

For the variance of mi, given i such that 2g ≥ i > 2g+1, we get a somewhat more complex expression

$${\sigma }_{i}^{2}={{{{{{{\rm{Var}}}}}}}}[{m}_{i}]=\frac{{\sigma }^{2}}{{4}^{g}}+\frac{({2}^{g}-1)\mu }{{4}^{g}}$$
(9)

For μi and σi, respectively, in (8) and (9), since i is constrained by 2g ≥ i > 2g+1, we can effectively substitute the cell label subscript i with the cell generation subscript g, such that for any integer value of i, the inequality 2g ≥ i > 2g+1 is satisfied

$${\mu }_{g}={{{{{{{\rm{E}}}}}}}}[{m}_{g}]=\frac{\mu }{{2}^{g}},$$
(10)
$${\sigma }_{g}^{2}={{{{{{{\rm{Var}}}}}}}}[{m}_{g}]=\frac{{\sigma }^{2}}{{4}^{g}}+\frac{({2}^{g}-1)\mu }{{4}^{g}}$$
(11)

Asymmetric case

If we assume \(p\,\ne\, \frac{1}{2}\), P(m2i∣mi) and P(m2i+1∣mi) have averages

$${\mu }_{2i}={{{{{{{\rm{E}}}}}}}}[{m}_{2i}| {m}_{i}]={m}_{i}p,$$
(12)

and

$${\mu }_{2i+1}={{{{{{{\rm{E}}}}}}}}[{m}_{2i+1}| {m}_{i}]={m}_{i}(1-p).$$
(13)

Similarly, for the variance, one finds

$${\sigma }_{2i}^{2}=E[{m}_{2i}^{2}]-E{[{m}_{2i}]}^{2}=p\cdot q\cdot {m}_{i}+{p}^{2}{\sigma }_{i}^{2}$$
(14)

and

$${\sigma }_{2i+1}^{2}=q\cdot p\cdot {m}_{i}+{q}^{2}{\sigma }_{i}^{2}$$
(15)

where q = 1 − p.

As the cell population grows, in an asymmetric division scenario, the different branches of daughter cells form subpopulations with distinct statistical expressions (see Fig. 3a). From these expressions, we understand how the mean and the variance of each branch in the lineage tree are linked to the mean and variance of the distribution of the initial progenitor cells population. In contrast to the symmetric case, in which the two branches of daughter cells originated by the same sub-population of mother cells form identical distributions, in the asymmetrical case, both the means and the variances are different. In terms of generations, starting with an initial population of progenitor cells at generation 0 with a localized distribution, we end up with two distributions at generation one, four at g = 2, and 2g at generation g (see Fig. 3b, c). In general, our capability to distinguish the different subpopulations of each generation depends on the distribution parameters. In our case, the branch sub-population variance is too large to distinguish the single population branches in the scatter plot Fig. 3b). Thus, to compare the model with the experimental data, for which we can only distinguish the distribution for each generation, we need to compute the mean and the variance of the cells at a certain generation, which is given by the superposition of all the lineage branches in the generation. Consequently, we want to compute for each generation

$${\mu }_{g}=\frac{1}{{2}^{g}}\mathop{\sum }\limits_{k=1}^{{2}^{g}}{\mu }_{g}^{k}$$
(16)

where \({\mu }_{g}^{k}\) is the mean of the k-th sub-population of generation g and we assume that in each generation all the cells that replicate give rise to two daughters.

Knowing the initial population mean value, μ, equation (16) can be solved (see Supplementary Note 1 for more details) as

$${\mu }_{g}={\left(\frac{1}{2}\right)}^{g}\mu$$
(17)

It is interesting to note that the dependence of p (the coefficient of the binomial partitioning) cancels out because we are considering the mean of the mean values.

Consequently, to estimate p we must compute the variance. In the case of the variance, we want the variance of a distribution that is given by the mixture of the distributions of all the daughter cells subpopulations. It can be demonstrated that the total variance is given by

$${\sigma }_{g}^{2}=\frac{1}{{2}^{g}}\mathop{\sum }\limits_{k=1}^{{2}^{g}}\left({\sigma }_{g,k}^{2}+{\mu }_{g,k}^{2}\right)-{\mu }_{g}^{2}$$
(18)

which after some computations, reported in Supplementary Note 1, yields

$${\sigma }_{g}^{2}=\alpha {\left(\frac{1}{2}\right)}^{g-1}\frac{1-{(2\beta )}^{g}}{1-2\beta }+{\beta }^{g}({\sigma }^{2}+{\mu }^{2})-{\left(\frac{1}{2}\right)}^{2g}{\mu }^{2}$$
(19)

with α = pqμ and \(\beta =\frac{{p}^{2}+{q}^{2}}{2}\). Figure 3d, e shows the behaviors of the mean fluorescence and its fluctuations as a function of the generations for different values of p. One can see that the trends are in line with those reported in Fig. 2c, d. In the next section, we show that the proposed model quantitatively explains the observations.

Comparison between model and data

We performed both a set of experiments with two-color stainings at a time (CTV-MITO and CTV-PKH) and another with three-color stainings (CTV-MITO-PKH). Working with two markers at a time ensured that during the sorting procedure a sufficiently large sample was collected. In fact, as for single marker flow cytometry to ensure peak detection, we sorted cells with a narrow gate centered at the dye intensity distribution mode. Consequently, the gate sorted a fraction f of the total population. In a multicolor experiment, if we assume that each marker is gated so that a fraction f of cells is sorted for each color, in total, we sort a fraction fd where d is the number of markers used. Thus, to achieve the same starting sample of a two-color experiment in a three-color experiment we had to prepare to sort a much larger population of cells than for a two-color experiment. Furthermore, with two markers, we set the collection gate and sorted all the stained cells in one sorting procedure (see Fig. 1b). While for a three-marker experiment we had to implement a subsequent sorting protocol with gates defined on the third-marker. Moreover, the two-color analysis allowed us to verify the consistency of the data because the CTV marker was present both in the two-color and three-color experiments.

In the Gaussian mixture clustering procedure, we analyzed two colors at the time, Fig. 4a shows how the mean changes when clustering was made on CTV vs MITO data, while Fig. 4b shows the mean values for each generation, considering CTV and PKH intensities together. We note that according to the binomial partitioning model the mean value of the marker intensities at each generation is expected to depend only on the mean of the initial population, μ, and on the generation, g (see Eq. (17)). Moreover, the μ dependence on g takes the form of a multiplicative constant, so that if we look at the mean intensity divided by the initial intensity as a function of the generation, we have a parameter-free way to confirm that the binomial model assumptions are at least partially correct. The CTV behaved identically in the two cases. Experimental data points are displayed with different markers, one kind for each experiment. Meanwhile, lines show the trends of the fluorescence mean values given by Eq. (17). Indeed, the intensity diminished by a factor of two passing from one generation to the next for all four experiments, in perfect agreement with the binomial predictions.

Fig. 4: Comparison between model and experiments.
figure 4

a, b From left to right, the mean intensity of the three markers, MITO (MitoTracker, mitochondria marker), CTV (CellTrace Violet, cytoplasm marker), and PKH (PKH Linker, membrane lipid marker) as a function of the generation, g, normalized over the mean intensity of the initial population. Different markers represent different experiments: stars, upper triangles, down triangles, and circles correspond to experiments one to four, respectively. Lines are the corresponding solutions of Eq. (17). CTV data has been analyzed twice with the Gaussian mixture method, one time associated with MITO and another with PKH. c, d Same as in (a, b) but with the variance of the intensity of each marker. Theoretical predictions (lines) are obtained from Eq. (19). For all experiments, the markers correspond to the mean of the values for each generation, while the sizes of the markers are comparable with the standard errors of the mean. e Pearson correlation between the intensities of MITO and CTV markers as a function of the generations. The gray dotted line represents the predictions of an uncorrelated binomial partitioning, starting with a correlated initial population. f Same as in (c) but with the correlation between PKH and CTV markers.

Full size image

However, considering the mean intensity of each generation does not allow us to evaluate if the partitioning process is symmetrical or not. In fact, as we can see in Eq. (17) and Fig. 3d, the resulting mean intensity is unaffected by the partition fraction, p. We note that if we were able to distinguish the single subpopulations, forming the distribution of each generation (see Fig. 3c), then we could study only the mean to fully characterize the partition mechanism. This can be seen in Eq. (12).

Since the total mean alone does not provide other information about the partition mechanism, we moved to consider the variance. The variance analytical expression is given in Eq. (19), and it explicitly depends on p. In particular, Fig. 4c, d displays the square root of the variance of the fluorescence divided by the square root of the variance at generation 0. Once more, data points are plotted as a function of the generation, and with different markers for each one of the four reported experiments. Furthermore, lines in the same figure represent the results of the best fits using the theoretical expression for the variance of Eq. (19). The expression depends on three parameters, i.e., the mean intensity at generation zero, μ, the initial variance, σ2, and the partition probability, p. Fixing the μs to those found in the analysis of the mean intensities, we are left with two free parameters. The partition probability, p, is of particular interest since it gives a measure of the asymmetry. Clearly, with p = 0.5 we get a completely symmetric partitioning process, and for p → 0, we get an asymmetric or biased partitioning process.

As for the mean, the variances were also found in excellent agreement with the predictions of the asymmetric binomial partitioning model for all the three cellular constituents we analyzed, although with different partition probabilities. In particular, cytoplasm appeared to undergo a symmetrical partition, with a p of ~0.50 both when measured in couple with mitochondria (Fig. 4c) and membrane proteins (Fig. 4d). Since the amount of cytoplasm is expected to correlate with the cell volume, its symmetrical partitioning is informative on the population size distribution and in accordance with what is known on Jurkat cell division. On the other hand, mitochondria and membrane proteins were both unevenly partitioned, i.e., they have a p significantly lower than 0.5. In particular, mitochondria fluorescence variance is well reproduced with p ~ 0.43, as shown in Fig. 4c. Similarly, a partitioning probability of ~0.36 gives the best description of the fluorescence variance associated with the PKH dye (Fig. 4d).

To check whether the results could be due to the particular choices of the fluorescent dyes and thus to avoid artifacts due to the different fluorescent photostabilities and properties, we performed additional experiments using different dye combinations to mark the cell membrane, cytoplasm, and mitochondria (see “Methods”). As one can see from Supplementary Figs. 7 and 8, results are robust to changes in the fluorescence dyes, even in the cases where there is fluorescence dye loss (see Supplementary Note 3, and Supplementary Fig. 3).

In addition, we want to stress that we were able to characterize mitochondria and membrane partitioning because they were associated with the cytoplasmic symmetrical marker. In fact, our observations show that symmetrical inherited dyes also have an additional practical interest because they can be used to support peak detection of noisier markers43. Indeed, while CTV is inherited in a highly symmetrical fashion in Jurkat cells and showed well-resolved division peaks, both the other dyes tended to partition asymmetrically. This effect broadens the fluorescence distributions. As one can see in the examples from Fig. 1e, g, where the single marker marginal distributions and the two-marker distributions are shown together, the two generations are hardly distinguishable if we look at the PKH intensity alone, while they are well separated if we use the CTV marker. Thus, the strategy of using different dyes enabled us to reduce the variation in the resolution peak for each generation and also to analyze correlations among different markers. In our case, the simultaneous staining of mitochondria and cytoplasm, and that of membrane and cytoplasm, enabled us to investigate the correlations between their partitioning mechanisms.

Specifically, we measured the Pearson correlation coefficient between the intensity of MITO and CTV and between PKH and CTV intensities in different generations. As one can see from Fig. 4e, f, the fluorescences of both couples tended to present an initial correlation higher than zero. The correlation of mitochondria and cytoplasm (MITO-CTV) is ≃0.35, while that of membrane proteins and cytoplasm is ρ ~ 0.20. During proliferation, both correlations were preserved from one generation to the next for all experiments where the initial distributions were correlated. Interestingly, the only experiment that started with an uncorrelated population (experiment 3 in Fig. 4f) showed a progressive increase of correlation which reaches the final value of the other experiments. Notably, the observed behavior of correlations is not compatible with uncorrelated two-variable binomial statistics. In fact, in this case, the expectation is a correlation that decreases over the generations even if the population of cells starts with a correlated distribution of stained compounds (the gray dotted lines in panels e and f of Fig. 4).

Discussion

Phenotypic heterogeneity has gained much attention in recent years as novel high-throughput experiments became able to study cell populations in greater detail. In particular, variability of non-genetic origin has been found to play a dominant contribution to cancer heterogeneity and to act as a driving force for tumor evolution47,48. While such variability is found in many biological aspects, its origin is often ascribed to stochasticity at the level of gene expression and its regulation49,50. However, fluctuations also arise during cell division when molecules are partitioned stochastically between the two daughter cells22,39,51,52. Disentangling the two contributions is quite difficult, and it has been suggested that often experimental outcomes may have been misinterpreted53. In fact, experiments often look at phenotypic traits to infer the effects of partitioning errors, even though those indirect observables are influenced by homeostatic pressure and thus by regulation. Furthermore, it is not clear how a homogeneous cell population is preserved despite partitioning errors or asymmetry. To systematically analyze these problems, we devised an experimental and computational apparatus to follow the proliferation of cell populations and directly measure how different cellular elements partition. Our method is based on staining cellular components with fluorescent markers at the beginning, sorting the stained cell population, and then measuring the dynamics of the partitioning process. This enables us to limit expression effects.

In general, the partitioning of cellular compounds is expected to depend on the particular types of compounds. Intuitively, those compounds that are present in huge numbers and are uniformly distributed across the whole cell could follow a random segregation mechanism. Meanwhile, both the replication and partitioning of nuclear chromosomes must be stringently controlled, with each chromosome being replicated exactly once per cell cycle and one copy transmitted to each daughter cell. To ensure perfect hereditary continuity of all the genes, stringent control is necessary to avoid the loss of essential genes in some cells. In fact, errors in the duplication and partitioning of genes or entire chromosomes are the cause of important genetic disorders, such as trisomy 2154. The scenario is less clear for cell organelles and internal compounds. Many RNAs, proteins, and organelles are present in numbers low enough that segregation of individual copies causes large fluctuations at cell division. Moreover, the presence of biases in the partitioning of some compounds can produce a systematic difference in the daughters, which fosters cell differentiation55. Even symmetrically dividing cells can then produce daughters with very different compositions by chance. Depending on the rates with which compounds are restored to their physiological levels, the single cell can generate variability. Most interestingly, as much as partitioning noise appears to take part in several cell division processes, cell populations are characterized by homeostasis that tends to contrast with noise56,57. Consequently, partitioning noise and other stochastic factors must cooperate to coordinate the population in a collective homeostatic state, similarly to other stochastic systems that present emergent collective properties58,59,60,61,62.

Multiple studies adopted the symmetric binomial partition model to describe partitioning because it naturally applies binomial distributions and requires the minimal assumptions that sub-components partition independently, and are inherited by daughter cells with equal probability54,63,64,65. This approach imposes that sub-components are segregated by a symmetrical process (symmetric segregation), and that uneven portions are the result of errors. The asymmetric partitioning model drops the equal probability assumptions and allows cellular sub-components to have a bias toward one of the daughter cells, with a single parameter p that measures the bias of the process and the noise and can be extrapolated from the trend of variance of the fluorescence distribution as a function of the generations. For p = 0.5 we have symmetric segregation and for p ≠ 0.5 we have asymmetric segregation. A recent work of Lijster and Åberg66 discusses the fact that to estimate the asymmetry parameter p from the fluorescence (or particle count) distribution we can not use the mean but propose to use the coefficient of variation of the full population. Unfortunately, the coefficient of variation is dependent on the relative frequencies of the generations of the cells in the population, and for this reason, it changes in time dynamically. With our analytical discussion, we show that under the asymmetric partitioning model assumptions, the variance of any single generation is constant. Furthermore, if we compute the variance of the cells in each generation, we can measure the asymmetry without taking into account the time of the measurements.

In particular, in the present work, we investigate the partitioning of three important cellular components during the proliferation of cancer Jurkat T-cell populations. We use three distinct fluorescent markers to label mitochondria, membrane proteins, and cytoplasm, and quantify both the degree of segregation asymmetry of each element between the two daughter cells, and the correlations in the partitioning of different kinds of compounds. Indeed, comparing the fluorescence mean intensities and their fluctuations against the predictions of the asymmetric binomial model, we found that mitochondria and double layer lipid membranes partition asymmetrically while cytoplasm partitions symmetrically. The degree of asymmetricity in mitochondria partition corresponds to a p ~ 0.43. It is important to consider that in comparison, to mammalian epithelial stem-like immortalized cells that asymmetrically partition more than five times the old mitochondria in the mother cell, for a Jurkat cell the bias is much smaller. This indicates that in Jurkat cells asymmetrical segregation is present even if in a weaker form in comparison to what is observed for stem-like asymmetrically dividing cells. Second, it shows how our method takes advantage of the large number of cells that can be processed in a flow cytometry experiment to improve the experimental resolution.

Similarly, a partitioning probability of ~0.36 gives the best description of the fluorescence variance associated with the PKH dye. In this case, the biased segregation may be linked to aggregates of membrane proteins (rafts) that T cells form when activated as well as from endocytosed fluorescence. Indeed, Jurkat cells are known to form clusters of membrane proteins. Those microdomains are required for protein-protein interaction and are implicated in T cells signaling during activation67. To verify that the membrane signal was due in some cases to clusters of aggregation, we acquired confocal microscopy images at the starting point (see Supplementary Fig. 9). Indeed, clusters of signals at the membrane level for the PKH fluorescence were evident in many of the analyzed cells.

Overall, our results have both technical and biological implications. From a methodological point of view, our method allowed us to overcome several technical issues. First of all, sorting the cells after the staining allows us to select a narrow initial distribution which is both essential to distinguish the generation peaks and very efficient for synchronizing the population. Second, Begum et al.43 show that succinimidyl fluorescent dyes can be used to track cell proliferation for a maximum of 72 h, while Filby et al.32 demonstrates that lipophilic dyes peaks are already hard to distinguish at 48 h, concluding that lipophilic dyes are weak proliferation trackers. Our multicolor approach allows us to resolve the peaks formed by the different generations up to 72 h for the lipophilic dye PKH, as for the CTV. This is due to the fact that with the use of multicolor flow cytometry, we were able to merge the information gathered from CTV on the analysis of the PKH. We can measure and analyze partitioning processes across different cell generations, and analyze how the covariance of the partitioning compounds changes. This allows us to show the gradual weakening of the PKH lipophilic dye peak resolution is the result of the partitioning noise that broadens the peaks. We also used the MitoTracker dye together with CTV and with both CTV and PKH, to track how progenitor mitochondria are inherited by the progeny. In conclusion, in contrast to a single-color flow cytometry live cell experiment that enables tracking cell proliferation of only one dye at a time32, in our multicolor experiment, we were able to simultaneously follow more than one live cell marker and quantify their partitioning during cellular proliferation up to more generations than one can do using only one marker at a time.

From a biological point of view, the high resolution of our method allowed us to observe asymmetric partition in Jurkat T-cell proliferation. It is worth noting, that even if the asymmetric division was observed in differentiating cells from cancer cell lines and in non-differentiating yeast cells, homogeneous cancer cell lines such as HeLa (similar to Jurkat) cells were used as an example of symmetrically dividing cells12. Being able to show that Jurkat cells preserve a small degree of asymmetric partitioning may have fundamental consequences on the way we conceptualize symmetric and asymmetric division because it implies that in non-differentiating cell lines asymmetric division persists (even if with a small value p). Furthermore, our present results can be more directly interpreted in the light of the inheritance of cellular sub-components that were part of the starting population stained at time zero. Indeed, the inheritance of marked organelles, indicates that these organelles were already present in the initial population, and are for this reason older organelles. Consequently, asymmetric division allows the same cells to inherit a larger portion of older organelles, while other cells inherit fewer organelles and are for this reason forced to produce newer organelles before they divide. While the biological role of asymmetric partitioning needs to be further investigated, we speculate that it could be connected to the anti-aging defense of the population as similar mechanisms have been observed in yeast cells and mammalian epithelial stem-like immortalized cells 4,18,19,21,68.

Finally, we want to stress that our apparatus can be easily applied to different cell lines and different kinds of compounds providing a powerful tool for understanding partitioning-driven heterogeneity. Moreover, our method is not based on transfection, and therefore in the future, it may be used in primary cell cultures and in in vivo studies. This is relevant since primary culture and in vivo studies are considered to be promising improvements on the traditional in vitro systems69. For instance, in future work, we could evaluate whether there exist differences between the partitioning of Jurkat T cells and of non-transformed T cells, for which the mechanisms that regulate growth are unaltered70. These kinds of comparisons could help us understand whether differences in partitioning play a role in tumor evolution. In perspective, this could be particularly relevant in the case of tumor microenvironment diversity, where comprehension of the high cell heterogeneity could pave the way to novel and more targeted therapies.

Methods

Cell culture

E6.1 Jurkat cells (kindly provided by Dr. Nadia Peragine, Department of Cellular Biotechnologies and Hematology, ‘Sapienza’ University) were used as a cell model for proliferation study and maintained in RPMI-1640 complete culture media containing 10% FBS, penicillin/streptomycin plus glutamine at 37 °C in 5% CO2. Cells were passaged three times prior to amplification for the experiment. Cells were then harvested, counted, and washed twice in serum-free solutions and resuspended in pre-warmed PBS for further staining.

Cells fluorescent dye labeling

To track cell proliferation by dye dilution establishing the daughter progeny of a completed cell cycle, cells were stained with PKH26, CellTrace™ Violet stain (CTV), and MitoTracker® Green FM according to the manufacturer’s instruction. The use of three different cell markers, respectively cell membrane, cytoplasm, and mitochondria, allowed for a multi-variable proliferation analysis. Fixable Viability Stain 780 (FVS780) was used in each experiment to determine cell viability and verify low dye toxicity. The first staining was performed with PKH26, a red fluorescent dye used for Cell Tracking Lipophilic Membrane Dyes (MINI26-1KT, Sigma, St. Louis, USA). Briefly, 20 × 106 cells were collected, washed once in PBS, and resuspended in 500 μl of Diluent C for lipophilic dyes. The cell suspension was then added to the same volume of 2x Dye Solution (4 μM) in Diluent C by adding 4 μl of the PKH26, mixed gently, and kept for 2 min in gentle agitation using vortex at room temperature (RT). Immediately after, the cellular suspension was blocked with an equal volume of complete FBS, added additional complete media to fill the tube, and centrifuged. Cells were washed again once in PBS prior to the following staining. The CTV dye (C34557, Life Technologies, Paisley, UK), typically used to monitor multiple cell generations, was performed according to the manufacturer’s instruction, diluting the CTV 1/1000 in 1 ml of PBS for 20 min at RT. Afterward, complete media was added to the cell suspension for additional 5 min incubation before the final washing in PBS. MitoTracker® Green FM (M7514, Molecular Probes, Eugene, USA), green-fluorescent mitochondrial stain, was used at the concentration of 20 nM (stock 20 μm in DMSO) in PBS for 30 min in a water bath at 37 °C, mixing every 10 min. Cells were then washed once in PBS. The live/dead staining, Fixable Viability Stain 780 (565388, BD Biosciences, San Jose, USA) was usually performed last, 1/1000 in PBS for 10 min at RT. For each experiment, a new aliquot of Viability marker was used. To confirm the data and to avoid artifacts due to the different fluorescent photostabilities and properties, as a control to the methodology, experiments were reperformed with different dye combinations to mark the cell membrane, cytoplasm, and mitochondria. For this purpose, in the following control experiments, we used the alternative dyes: CellTrace™ Yellow (C34573 A, Life Technologies, Paisley, UK), CellTrace™ Far Red (C34572, Life Technologies, Paisley, UK), MitoTracker® Deep Red 633 FM (M22426, Molecular Probes, Eugene, USA) and PKH67 green-fluorescent cell linker (MINI67-1KT, Merck Life Science, Milano, Italy).

Cell sorting

Jurkat cells labeled with dyes were sorted using a FACSAriaIII (Becton Dickinson, BD Biosciences, USA) equipped with Near UV 375, 488, 561, and 633 nm lasers and FACSDiva software (BD Biosciences version 6.1.3). Data were analyzed using FlowJo software (Tree Star, version 9.3.2 and 10.7.1). Briefly, cells were first gated on live cells, viability dye exclusion, followed by doublets exclusion with morphology parameters, both side and forward scatter, area versus width (A versus W). The unstained sample was used to set the background fluorescence for each channel. For each fluorochrome, a sorting gate was set around the max peak of fluorescence of the dye distribution32. In this way, the collected cells were enriched for the highest fluorescence intensity for all the markers used according to the experimental setup: PKH, CellTrace, and/or MitoTracker. Following isolation, an aliquot of the sorted cells was analyzed at the same instrument to determine the post sorting purity and population width, resulting in an enrichment >99% for each sample (see Supplementary Fig. 4).

Cell culture for dye dilution assessment

The sorted cell population was seeded in 6 wells plates (BD Falcon) at 1 × 106 cells/well and kept in culture for up to 72 h. To monitor multiple cell division, an aliquot of the cells in culture was analyzed every 18, 24, 36, 48, 60, and 72 h for the expression of the different markers by the LSRFortessa flow cytometer. In order to set the time zero of the kinetic, prior culturing, a tiny aliquot of the collected cells was analyzed immediately after sorting at the flow cytometer. The unstained sample was used to set the background fluorescence as described above. Every time that an aliquot of cells was collected for analysis, the same volume of fresh media was replaced with culture.

Expectation–maximization and the Gaussian mixture model

We used the expectation–maximization (EM) algorithm to detect the clusters in Gaussian mixture models71. The EM algorithm is composed of two steps: the expectation (E) step and the maximization (M) step. In the E-step, for each data point f, we used our current guess of πg, μg, and σg, to estimate the posterior probability that each cell belongs to generation g given that its fluorescence intensity measure as f, γg = P(g∣f). In the M-step, we use the fact that the gradient of the log-likelihood of p(fi) for πg, μg, and σg can be computed. Consequently, the expression of the optimal value of πg, μg, and σg is dependent on γg. It is shown that under, certain smoothness condition the iterative computation of E-stem and M-step leads us to the locally optimal estimate of the parameters πg, μg, and σg, and returns the posterior probability γg which weights how much each point belongs to one of the clusters. Here, we used this model to perform cluster analysis and detect the peaks which correspond to different generations. Then, we estimated πg, E[fg], and Var[fg] from these clusters.

Gaussian mixture clustering

To estimate to which generation the detected cell belongs, we used the Gaussian mixture method. For the experiments with two dyes (CTV-PKH and CTV-MITO) we always used all two fluorescence measures for the clustering method. In the three dyes case (CTV, PKH, and MITO), we used all three fluorescence measures, except when the MITO staining was getting close to the background fluorescence in the successive generation, since this was biasing the results of the clustering method. Figure 1d and Supplementary Fig. 5 show representative examples of how the Gaussian mixture method was able to cluster the data. For experiment 1, we recorded 3546 cells with the cell sorter at time 0, and 10,000 cells with the flow cytometer at times 24, 36, 48, 60, and 72. For experiment 2, we recorded 4811 cells with the cell sorter at time 0, and 10,000 cells with the flow cytometer at times 18, 24, 36, 48, 60, and 72. For experiment 3, we recorded 1537 cells with the cell sorter at time 0, and 10,000 cells with the flow cytometer at times 18, 24, 36, 48, 60, and 72. For experiment 4, we recorded 1000 cells with the cell sorter at times 0, 18, 24, 36, 48, 60, and 72. Outlier cells that biased the estimation of the generation peaks were dropped.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.