Mitochondrial heterogeneity, metabolic scaling and cell death

Heterogeneity in mitochondrial content has been previously suggested as a major contributor to cellular noise, with multiple studies indicating its direct involvement in biomedically important cellular phenomena. A recently published dataset explored the connection between mitochondrial functionality and cell physiology, where a non‐linearity between mitochondrial functionality and cell size was found. Using mathematical models, we suggest that a combination of metabolic scaling and a simple model of cell death may account for these observations. However, our findings also suggest the existence of alternative competing hypotheses, such as a non‐linearity between cell death and cell size. While we find that the proposed non‐linear coupling between mitochondrial functionality and cell size provides a compelling alternative to previous attempts to link mitochondrial heterogeneity and cell physiology, we emphasise the need to account for alternative causal variables, including cell cycle, size, mitochondrial density and death, in future studies of mitochondrial physiology.


Introduction
Non-genetic heterogeneity describes observable differences between genetically identical cells, and is a fundamental observation in cell biology [1][2][3], with implications for important biological phenomena including development [4,5], bacterial persistence [6], cancer malignancy and drug resistance [7]. It is becoming increasingly recognised that variability in cellular populations of mitochondria -organelles of central bioenergetic importance -is a significant source of cellular heterogeneity in eukaryotes. As indicated experimentally by Das Neves et al. [8], and theoretically by Johnston et al. [9], mitochondrial variability could be a major driver of noise in global transcription rate which, in turn, is a key contributor to the extrinsic cellular noise observed in a wide range of downstream processes [10][11][12][13][14]. Indeed, recent studies have suggested that variability in mitochondrial content plays a direct role in various stochastic cellular processes, such as differentiation [15,16], cancer metastatic potential [17] and chemotherapeutic resistance [18]; metabolic heterogeneity has also been observed to underlie fluctuations in prokaryotic systems [19]. These associations underscore the importance of continuing to expand our understanding of the origins and consequences of mitochondrial heterogeneity.
The sources of mitochondrial variability are diverse. For instance, mitochondrial membrane potential (DC) is heterogeneous in both intracellular mitochondrial populations [20] and at the aggregate level between cells [8]. For intracellular mitochondrial populations, organelles with lower DC are less likely to join the mitochondrial network and tend to be eliminated by mitophagy [20], suggesting that DC is a measure of mitochondrial functionality. Disruption of mitochondrial fusion can also lead to intracellular heterogeneity in DC [21]. While the precise function of mitochondrial networks is unclear [22], it is likely that a wide range of cellular functions will depend non-linearly on the degree of mitochondrial fusion. It is therefore reasonable to speculate that cell-to-cell variability in the mitochondrial network state may have non-trivial consequences for cellular heterogeneity, for instance, via intermediates such as ATP or reactive oxygen species production.
Since mitochondria possess their own genomes, which are replicated, degraded and stochastically partitioned during cell division [23], variability in copy number and sequence of mitochondrial DNA (mtDNA) are also important sources of bioenergetic variability, most dramatically illustrated by the profound phenotypes associated with mtDNA diseases [24]. Given the existence of two co-existing mtDNA variants, and stable mean copy numbers of each species, mathematical modelling suggests that the variance in mutant load (or heteroplasmy) of mitochondrial DNA inevitably increases over time [25]. Thus, given the existence of more than one mtDNA variant, mtDNA heteroplasmy is also a contributor to mitochondrial variability.
Recent work by Miettinen and Bj€ orklund [26] has generated a rich dataset probing the relationship between mitochondrial functionality and cell size, as an additional axis of cellular heterogeneity. The authors also explore the connection with other cellular features such as apoptosis, cell cycle and metabolism. In this paper, we consider a set of quantitative arguments that could account for this dataset. We invoke several models for bioenergetic scaling, the dynamics of cell proliferation and death, and the coupling between cell physiology and mitochondrial content. Amongst the insights gained from these models, we find that: Cell size scaling of cellular power demand may explain the observed reduction in mitochondrial functionality with increasing cell size, in large cells; The cell's characteristic (or 'optimal') size may be largely set by a particular form of intrinsic cell growth dynamics, rather than mitochondrial functionality; Cell death alone may be used to explain the observed nonlinearity in mitochondrial functionality with cell size, if mitochondrial functionality is considered a passive indicator of the health of a cell; A non-heritable, non-linear relationship between mitochondrial functionality and cell size, as illustrated by Miettinen and Bj€ orklund, is compatible with a wider set of single-cell data [8] and provides a compelling alternative to link mitochondrial heterogeneity and cell physiology.
By demonstrating the existence of multiple competing hypotheses which may account for the data of Miettinen and Bj€ orklund [26], we highlight the importance of considering variables such as cell cycle, volume, mitochondrial density and cell death in future studies of mitochondrial heterogeneity. We propose further experiments to distinguish these models and further elucidate the coupling between mitochondria and cell dynamics.

Linking mitochondrial heterogeneity and cell size
Mitochondrial membrane potential is often considered a measure of mitochondrial functionality, since a higher membrane potential may support a larger respiratory rate and ATP/ADP ratio [27,28], as well as its role in quality control [20] (although this correspondence may not hold under all conditions, since inhibiting ATP synthase can cause membrane potential to increase [28,29]). It is important to distinguish between total DC, which is expected to scale with a characteristic size of the system (e.g. the size of the mitochondrial population, total cellular protein or some measure of the size of the cell itself, e.g. the radius), and a measure of DC which is normalised to account for this scaling. We denote a relative measure of DC with rDC m , rDC p or rDC r for normalisations by mitochondrial mass, total cellular protein and cell radius, respectively.
The main probe which Miettinen and Bj€ orklund [26] used for normalising DC was the lipophilic cationic dye JC-1 [30]. JC-1 can pass cellular membranes and accumulates in regions of negative potential, such as the mitochondrial matrix. JC-1 is green-emitting in its monomeric form; however, if its concentration exceeds a threshold value then it forms aggregates which are red-emitting [30]. By taking the ratio of red to green signals, JC-1 is considered as a semi-quantitative measure of rDC m . Perry et al. [29] have suggested several complexities in using this probe. Since aggregates form according to threshold effects based on the concentration of monomers, JC-1 may be ill-suited to detecting subtle gradations in DC and has been recommended to be used as a binary indicator of high/low DC. Furthermore, it has been shown that aggregate formation is sensitive to loading concentration [31], suggesting that cellular surface-to-volume ratio may have an effect on aggregation rate. This is of particular relevance to the data of Miettinen and Bj€ orklund [26] since variations in cell size are of interest. Finally, equilibration times tend to be long for JC-1, as can be seen in [26] Fig. S2E (although the qualitative observations of this particular experiment appeared to hold throughout the staining kinetics).
Despite these complexities, many of the key insights of Miettinen and Bj€ orklund [26] were reproduced with other means of normalising mitochondrial membrane potential. Throughout their experiments, the authors used flow cytometry to measure cell size using forward scatter (FSC). Calibration of FSC is important for its rigorous interpretation [32]. Upon binning cells by diameter, with bin width of $100 nm, the authors found good correlation between FSC and cell radius across the binned populations. DC was measured with TMRE (another lipophilic cationic dye) and total cellular protein content using CellTrace FarRed stain [33] which covalently binds to cytosolic proteins. The ratio of these signals yields rDC p . The authors found that rDC p initially increased and subsequently decreased (henceforth referred to as 'turning behaviour') with increasing cell radius, which was qualitatively similar to rDC m found via JC-1 ( [26] Fig. S2K, 1C). The authors also used MitoTracker Red as a measure of DC and MitoTracker Green as a measure of mitochondrial mass ( [26] Fig. 1D). Their data suggested that rDC r shows turning behaviour with increasing cell radius; reanalysis of these data also suggests that rDC m shows turning behaviour with either increasing cell radius or volume (see Supporting Information Fig. S1).
While the precise interpretation of JC-1 ratio remains somewhat unclear due to the complexities above, the empirical observation that median JC-1 monomer fluorescence closely tracks median total cell protein suggests that it is appropriate for normalisation by cell size [26]. The orthogonal validation experiments suggest some degree of qualitative similarity between rDC m , rDC p and rDC r , all of which show turning behaviour with increasing cell radius. The qualitative agreement between median JC-1 ratio and these validation experiments suggest that JC-1 ratio is a qualitative measure of system-size-normalised mitochondrial membrane potential. We, therefore, use the notation rDC to denote JC-1 measurements, which the data of [26] suggests has qualitative similarity to rDC m , rDC p and rDC r . This turning behaviour was found with the JC-1 probe in immortal cell lines from human (Jurkat), fly (Kc167) and chicken (DT40) as well as primary rat hepatocytes and human primary HUVECs (i.e. non-immortalized cells, [26] Fig. S2L), providing some evidence that this observation may hold beyond immortalized cell cultures.
Miettinen and Bj€ orklund [26] found that the maximum rDC over cell size approximately coincided with the modal cell size -see Fig. 1. The distribution of rDC was also invariant to cell cycle stage for cells of an approximately fixed size. The authors therefore postulated that non-linearity in mitochondrial functionality results in an 'optimal' cell size [26], independently of cell cycle stage.

Cell physiological models coupled to cell death
Miettinen and Bj€ orklund speculate that an allometric decline causes reduction in mitochondrial functionality with increasing cell size [26,34], citing the West, Brown and Enquist (WBE) model of allometric scaling of metabolism with body size [35]. Briefly, this model suggests that transport of materials through fractal space-filling networks explains the relationship between organismal mass and metabolic activity (see Box 1 for a discussion of this hypothesis). In contrast, the authors speculate that small cells display low rDC due to a 'newborn effect' whereby daughter cells inherit the low mitochondrial functionality of their parents, which itself originates from the allometric decline in the final stages of growth in parent cells. Intermediate-sized cells were postulated to 'reset' their metabolic activity, increasing it during early stages of growth [26,34]. However, both the 'newborn' and metabolic 'reset' arguments suggest differences in rDC at different stages of the cell cycle, which was not observed using JC-1 ( [26] Fig. 2D, S3C and D). We note that the lack of cell-cycle dependence in rDC was not reproduced independently of the qualitative JC-1 probe, so the newborn/reset hypothesis remains possible.
Here, we take the approach of first considering a combination of metabolic scaling and cell death as a possible link between mitochondrial behaviour and cell physiology. We will later consider an alternative purely involving cell death as the underlying causal variable of the observed non-linearity.
A bioenergetic proposal: A combination of scaling of cellular power demand with cell size, and cell death, may account for observed patterns of mitochondrial functionality In Box 2, we propose that metabolic scaling with a particular model of cell death may recover the turning behaviour in rDC with cell radius observed by Miettinen and Bj€ orklund [26]. The model suggests that, under non-pathophysiological circumstances, cells alter their mitochondrial functionality to maintain energy supply/demand balance at different cell volumes (v).
The model makes the non-intuitive assumption that cells do not alter their mitochondrial density in response to volumedependent changes in cellular energy demands: one may speculate that this is due to the timescale of changes in rDC being able to match cellular energy demands, whereas mitochondrial mass (n) may not. Proportionality between n and v has found some support in the literature [37][38][39]. Indeed, recent mathematical modelling [46] of the cellphysiological consequences of heteroplasmy in a deleterious For each bin, the authors measured total mitochondrial membrane potential, normalised by cell size, (rDC) using the ratiometric probe JC-1, and reported its median (red curve). A: Approximate coincidence between the mode of the cell radius distribution (grey histogram), determined through forward scatter with flow cytometry, and median rDC. B: Although cells were gated for cell death using propidium iodide, the smallest cells in the population were high in an alternative annexin-based sensor of cell death (pSIVA) after normalisation by cell protein content, suggesting that these were early apoptotic cells. Observations A and B do not exclusively rely on JC-1 measurements. C: For approximately fixed cell size, the authors showed that the distribution of rDC was approximately invariant across different cell cycle stages (red, blue, orange histograms), suggesting that the variation in rDC is more strongly informed by cell volume than cell cycle stage. Note that this observation was not validated independently of JC-1. Based on figures from [26]. mtDNA mutation associated with MELAS [47] has suggested that failure to maintain the cytoplasmic density of wild-type mitochondrial DNA is associated with a pathophysiological response in vitro. However, the data of [26] (neglecting JC-1 measurements) did not strongly constrain the precise relationship between n and v so other scaling relationships are feasible, see Supporting Information Section S2 for further discussion.
The scaling argument culminating in Equation (3) suggests that rDC decreases with cell volume. This is because mitochondrial mass has been taken to scale proportionately with volume but power demands do not (Equation (1)). Small cells have lower mitochondrial mass, and must therefore have higher mitochondrial functionality to satisfy the constant component of power demand placed by, e.g. the nucleus (k n ). As cells become larger and have higher mitochondrial mass, power demand from the nucleus may be more easily satisfied, so mitochondrial functionality reduces towards a constant value k v .
This model alone suggests that rDC is a purely decreasing function of cell volume: however, Miettinen and Bj€ orklund observed a turning point in the relationship between rDC and v. The data in [26] show increased cell death in the smallest cells in the population [26] (Fig. 3C), and potentially also a subtle increase in the largest cells. Furthermore, the authors observed a depolarized subpopulation of cells at low JC-1 monomer values ( [26] Fig. S2G), potentially indicative of a subset of cells with low mitochondrial functionality at low cell sizes. Cell shrinkage is associated with apoptosis [48], as is loss of mitochondrial membrane potential [49]. Hence, we suggest that smaller cells have a higher cell death rate, and that mitochondrial functionality is a passive indicator of cell death under these circumstances ( Fig. 2A). The result of simulating this model in a proliferating population of cells undergoing 'adder' dynamics [50,51], whereby cells add a constant amount of volume each cell cycle, is presented in Fig. 2.
As can be seen in the distribution of cell volumes (Fig. 2C), adder dynamics imposes a characteristic cell size associated with the mode of the distribution, which is often referred to as an 'optimal' cell size. The width of the final distribution is set by the variability in partitioning of volume between daughters at division (although other sources of variability may be envisaged such as noise in the growth rate or the amount of volume added per cell cycle [54], which are neglected here). Under these simulation dynamics, the modal cell volume and the maximum value of rDC are essentially independent (Fig. 2C) and can be tuned relatively to each other through appropriate choices of model parameters. The lack of alignment between the modal cell size and maximal rDC could also be achieved with 'sizer' dynamics, where cells divide once exceeding a maximal cell size. The coefficient of variation (CV) in rDC also showed a minimum (Fig. 2D), in qualitative agreement with [26], although the minimum CV was constrained to coincide with the minimum rDC (see Supporting Information Section S4 for further discussion). It is therefore largely compatible with the data of [26] to suggest that a cell's characteristic size is set by a particular form of intrinsic cell dynamics, rather than requiring a direct link with rDC as proposed in [26].
In Supporting Information Section S5, we explore an alternative hypothesis where cell death rate is a non-linear function of cell size, and that mitochondrial membrane potential is a passive readout of cell death rate over all values of cell volume. We found that this is also a compatible interpretation of the data presented in [26]. However, the causality of this model may also be inverted, where Box 1

Allometric scaling hypothesis of cellular power demands
The WBE model of allometric scaling [35] claims that the power-law relationship linking an organism's metabolic rate (B) and mass (M) is B / M 3/4 and can be derived by considering the transport of materials through spacefilling fractal networks of branching tubes, for instance the cardiovascular network of most vertebrates. The theory was subsequently extended to explain the metabolic requirements of single cells in terms of mitochondrial output [36]. Assuming that cell mass and volume (v) are proportional, through constant cellular density, and that the power demand (D) of a cell equates to its metabolic output, the WBE model suggests that the power demand of single cells scales as D / v 3/4 . We model mitochondrial power supply (S) in terms of mitochondrial mass density (r ¼ n/v, for mitochondrial mass n) and mitochondrial functionality (f) as S / rvf. Miettinen and Bj€ orklund [26], as well as other authors [37][38][39], suggest that n / v (i.e. r ¼ const) and therefore S / v(rDC) (see Main Text for further discussion of r variation with v), where we have interpreted f / rDC. Using S ¼ D results in the prediction that rDC / v À1/4 , in contrast to the model we suggest in Equation (3).
Over the relatively small dynamic range of cell volume in a population of cells, the rDC / v À1/4 relationship of the WBE model and Equation (3) are unlikely to make large differences in their predictions. However, there are some conceptual difficulties in accepting the WBE model above a simpler energetic scaling argument. Whilst mitochondria do form fractallike networks [40], and the level of fusion of these networks might have non-linear effects on metabolic output [22], it is not clear that a system of static, hierarchical, tubes with fluid flow is an appropriate modelling framework to describe energy transportation within single cells. Because of this, and the criticisms sometimes levelled at allometric scaling [41][42][43], we have considered a scaling argument which focusses on the potential bioenergetic demands from different specific cellular compartments [44]. Nevertheless, the fusion status of the mitochondrial network remains an appealing potential explanation for the mechanism by which cells alter their mitochondrial functionality to maintain energetic homeostasis with varying cell volume. It therefore remains an interesting area of further study to quantify the extent of mitochondrial fusion and fission [45] as a function of cell volume to explain the rDC variations observed by Miettinen and Bj€ orklund [26]. mitochondrial functionality is a non-linear function of cell volume, and cell death rate is linearly related to mitochondrial functionality. This interpretation aligns more closely with that of Miettinen and Bj€ orklund [26], who suggest that the non-linearity between rDC and v underlies their observations. This highlights the importance of excluding cell death as a potential confounding variable in future confirmatory studies so that these hypotheses may be disambiguated.
Future verification studies may include gating cells in FACS analysis based on the pSIVA cell death probe used in [26]. Alternatively, fluorescence labelling of upstream apoptosis markers, such as cytochrome-c, may allow detection of cells in early stages of cell death [55]. Combining this with fluorescence microscopy and tracking mitochondrial membrane potential, mitochondrial mass and cell volume, may allow early-apoptotic cells to be excluded and verification of the non-linearity between cell volume and relative mitochondrial membrane potential.
A parabolic function of mitochondrial functionality against cell size is compatible with a mathematical model of mitochondrial variability A previous model, accounting for the data in [8], described the cellular consequences of mitochondrial variability in HeLa cells [9]. That model phenomenologically described mitochondrial functionality as a stochastically inherited quantity, determined through a mean-reverting process (an AR(1) process). The biological mechanism for the heritability of mitochondrial functionality was unclear from that study, although one may speculate that the inherited concentration of mitochondrial fusion proteins may be a feasible explanation for this. Indeed, recent evidence suggests that the Power demand scaling (black dotted) was determined through a power supply ¼ demand relationship in Equation (3). Cell death scaling (black dashed) was determined through a piecewise linear relationship in volume (Equation (5)). The mean rDC (grey) was determined by taking the minimum of the demand scaling and cell death curves at every point in v (to model switching behaviour which is continuous in v). Addition of Gaussian noise to the mean relationship recovers a joint density of v and rDC (red shaded region) similar to the findings of Miettinen and Bj€ orklund ( [26] Fig. 1B). C and D: Cells were binned by volume in silico, and for each bin the mean and standard deviation of rDC and v were computed. C: Mean rDC (red points) shows a maximum with respect to volume which is not aligned with the modal cell volume (grey histogram of cell volumes), recapitulating the observation of Miettinen and Bj€ orklund ( [26] Fig. 1C). D: The coefficient of variation of rDC (red points) shows a minimum with respect to volume, as observed by Miettinen and Bj€ orklund ( [26] Fig. 3A), although the location of the minimum CV and maximum average rDC are constrained to coincide in our model (C and D) -see Supporting Information Section S4 for further discussion. For simulation details, see Supporting Information Section S3. proliferation rate of yeast is heritable, and that slow-growing cells have higher oxidative stress whose proliferation rate may be rescued with an anti-oxidant, suggesting a link to mitochondrial functionality [56]. Further studies are required to determine whether this is true for a larger class of systems, including mammalian cells.
The data of Miettinen and Bj€ orklund [26] immediately suggests an alternative nonlinear relationship between function and dynamics (see Fig. 1). We hypothesised that this relationship may be sufficient to account for the data in [8] without invoking a phenomenological mechanism for the stochastic inheritance of mitochondrial functionality (though the stochastic inheritance of mitochondrial mass remains supported). Briefly, the model of Johnston et al. [9] predicts the relationship between cell volume, mitochondrial content and mitochondrial functionality with a pair of coupled ordinary differential equations. Cells undergo sizer dynamics, and the population size is controlled to be static over time.
Here, we replace the stochastic inheritance of mitochondrial functionality by approximating the result of Miettinen Box 2

Mathematical model of power demand scaling and cell death
We will relate mitochondrial functionality (rDC(v)) to cellular volume (v) by modelling rDC as a combination of two models: the first describes the eventual decrease of rDC(v) with v using metabolic scaling arguments, denoted by rDC s . In contrast, the initial increasing part of rDC(v) will be modelled using cell death arguments and denoted by rDC m . A combination of rDC s and rDC m will be used to represent the turning behaviour observed in rDC with cell volume and therefore also cell radius (as the two are related through a monotonic relationship).
We relate mitochondrial functionality (rDC) to cellular volume (v) by considering total cellular power supply (S) and demand (D). To represent D, we consider contributions which scale with cell volume (e.g. cytoplasmic protein content), v, surface area (e.g. work done at the plasma membrane), v 2/3 , and constant terms which do not scale with cell size (e.g. replication of the genome). This may be written as where k i are positive real constants, which were loosely based on a study of the power demands of thymocytes [44] (see Supporting Information Section S3 for further discussion). In general, we might expect S to be of the form S / vr(v)f(v), where r(v) is mitochondrial mass density (r(v) ¼ n(v)/v, where n(v) is mitochondrial mass) and f(v) is a measure of mitochondrial functionality, in an analogous form to the bioenergetic driving term in [9]. We interpret f(v) / rDC(v), since larger values of rDC support a larger respiratory rate and ATP/ADP ratio [27,28]. We take r (v) ¼ const for parsimony (see Main Text for discussion). Hence, We model power supply as being balanced by power demand, S ¼ D, with cells varying rDC to maintain this relationship for various v. Rearrangement for rDC yields where the subscript s denotes this metabolic scaling model for rDC. For the cell death component of mitochondrial functionality, we assume that the cell death rate (m(v), the probability of cell death per unit time of individual cells) increases amongst cells with smaller volumes. Above a threshold volume (determined by continuity), cell death is modelled as negligible. We then model relative membrane potential as a passive indicator of the cellular death rate m, since cells which are apoptotic tend to have low mitochondrial membrane potential [49]. We hence use the following linear functions, where c m , c c , m m and m c are positive real constants. With this, we recover an alternative description of rDC which initially increases with volume because rDC m , m and v are each connected by negative linear relationships. We combine rDC s and rDC m through a rule which ensures continuity in rDC. For mathematical convenience, we choose rDC(v) ¼ min{rDC s (v), rDC m (v)} for all v. This function will show a turning point as cells switch from having their membrane potential governed by risk of cell death to power demand scaling, for increasing cell size. rDC is then observed with Gaussian noise N(0,s c ) with standard deviation (SD) s c , which accounts for technical and biological sources of noise.
We used a stochastic hybrid systems approach [52] to simulate a population of exponentially growing cells undergoing 'adder' dynamics [50,51], whereby cells add a constant amount of volume (K v ) to their initial volume (v 0 ) each cell cycle, which has found recent support in bacteria [53]. Cells are modelled to then divide by stochastically partitioning their volume, with the first daughter inheriting volume v 1 ¼ N((v 0 þ K v )/2,s v ) (with SD s v ), and the second inheriting v 2 ¼ v 0 þ K v À v 1 through conservation of mass. Cell death events occur stochastically according to the rate m(v) (Equation (4) and Bj€ orklund as a non-heritable parabolic relationship between mitochondrial functionality (which we interpret as rDC) and cell volume v where f max is the peak value of f, occurring when v ¼ v opt , and f m is the steepness of the parabola. In implementing the original model from [9], it was found that constraining f > 0.1 from the AR(1) process (rather than f > 0.01 described in [9]) yielded a better comparison to experimental cell cycle lengths. The result of replacing the AR(1) process in [9] with the deterministic relationship in Equation (6) is shown in Fig. 3. We found that whilst the distributions of cell volume and mitochondrial mass were relatively unchanged ( Fig. 3A and B), mitochondrial functionality had greater skewness under the parabolic model (Fig. 3C). The distribution of total mitochondrial content nf was narrower in the parabolic model (SD ¼ 11.5) than the original model (SD ¼ 17.4) (Fig. 3D) which has poorer agreement with data from [8]. However, the cell cycle length distribution in the parabolic model had much better agreement in its standard deviation (SD ¼ 5.07) than the original model (SD ¼ 25.2) compared with data (SD ¼ 5.13) (Fig. 3H). The low degree of correlation between cell cycle lengths between parent and daughter cells was also retained (Original: R 2 ¼ 0.11, Parabolic: R 2 ¼ 0.05) in agreement with the data of [8] ( Fig. 3C of [9]). Overall, this shows that the non-heritable non-linear mechanism between mitochondrial functionality and cell volume found in [26] may broadly account for the data in [8] without invoking a phenomenological stochastic model for the inheritance of mitochondrial functionality.

Conclusion
Heterogeneity in mitochondrial content is exerted along many axes, including intrinsic functional variability, network status and genetic variability. A host of cellular phenomena are affected by noise, a growing number of which are being directly attributed to mitochondrial variability. The study of Miettinen and Bj€ orklund [26] has extended this picture by probing the link between mitochondrial content and cell volume, suggesting a non-linear relationship between mitochondrial functionality and cell radius.
We have shown that a simple model of power demand scaling predicts a decreasing trend in relative mitochondrial functionality with increasing cell volume. By combining this with a simple cell death model, we were able to qualitatively account for a large number of observations by Miettinen and Bj€ orklund [26]. We also found that incorporating the nonlinear relationship between rDC and cell volume as suggested by Miettinen and Bj€ orklund into a model of mitochondrial heterogeneity [9] was broadly consonant with a wider set of single cell data [8]. We have found that it is also compatible with the data of [26] to suggest that cell death is the causal Figure 3. Nonlinear connection between mitochondrial functionality and cell volume is consonant with a wider set of single cell data. Replacing the stochastic process determining mitochondrial functionality (Original) with the parabolic relationship in Equation (6) (Parabolic), and comparisons with data from [8] where available. A: Distribution of cell volumes is relatively unchanged by the modification of f. B: Distribution of mitochondrial mass in both models are comparable with data. C: Distribution of mitochondrial functionality. Since Equation (6) is deterministic, f cannot exceed its maximal value f max , causing the skewness observed. D: The distribution of total mitochondrial functionality nf has a reduced variance in the parabolic model, which departs somewhat from the data of [8]. E: Relationship between ratio of cellular volumes at birth and the ratio of cell cycle lengths for sister pairs. F: Relationship between ratio of mitochondrial contents at birth and the ratio of cell cycle lengths for sister pairs. Both models are comparable with data in E and F. G: Mean cell cycle lengths AE standard deviation for the original (O) and parabolic (P) models, as well as the experimental observation in [8] (Data). The variance in cell cycle lengths is greatly reduced, since an additional source of variability is removed in P from the deterministic relationship Equation (6), which has better agreement with data. variable in determining the non-linearity in rDC. Our analysis has shown the existence of multiple competing hypotheses which are able to account for the data of Miettinen and Bj€ orklund, thus highlighting the importance of accounting for cell cycle, volume, mitochondrial density and cell death in future single-cell studies of mitochondrial heterogeneity.