Mitochondrial Variability as a Source of Extrinsic Cellular Noise

Model

While the heterogeneity of mitochondria has been observed experimentally and connected to variability in processes like transcription [31] and stem cell differentiation [38], the mechanisms by which mitochondrial variability influences other cellular processes has not been elucidated clearly. Here, we describe a simple model which formalises these links, and note that it is consistent with recent experimental results concerning mitochondrial heterogeneity (and variability in connected cellular features) [31]. The simplicity of our model means that analytic expressions can be derived for many quantities of interest, facilitating a more complete and intuitive understanding of the modelled biological connections. We will then use this model to investigate more specific questions regarding the links between mitochondrial variability and transcription rate and stem cell differentiation.

The central concept behind our model is illustrated in Fig. 1. Individual cells are characterised by three key variables: the volume of the cell (); the amount of mitochondrial mass in the cell (); and the degree of mitochondrial functionality (). This last quantity, , represents a coarse-grained measure of the efficiency of mitochondria within a cell – a factor which may be affected, for example, by the levels of reactive oxygen species (ROS), mitochondrial membrane potential, variability in mitochondrial protein complex abundance, and genetic differences between mitochondria [48]. ATP concentration in the cell is modelled as a function of these three quantities, and transcription rate is modelled as a function of ATP concentration [31]. Variability in cell volume, mitochondrial mass and mitochondrial functionality arises due to stochastic inheritance of these quantities at cell divisions. This variability causes cell-to-cell differences in ATP levels, and hence transcription rate, in a population of cells.

and transcription rate

We are interested in the time evolution of as a potential stochastic influence on downstream processes. Ref. [31] found ATP levels in the cell to be proportional to mitochondrial mass () and membrane potential (a factor that may be absorbed into our measure of ‘mitochondrial function’ ), motivating our choice of expression for ATP concentration:

(1)

In this expression, is a constant of proportionality linking the quantities within our model to a biological ATP concentration, and the meaning of the variable now becomes apparent as a scalar multiple relating mitochondrial density to . We note that other choices for the form of , including ODEs, are possible, and explore some alternatives in ‘Other Models’ (Text S1). das Neves et al. also show a link between the total transcription rate in a cell (measured through bromo-uridine incorporation across the whole nucleus) and , a sigmoidal curve, which we approximate (see ‘Parameterisation of ’ in Text S1) with

(2)

das Neves et al. record a change in the structure of this sigmoid curve in experiments where cellular chromatin is artifically decondensed. In these situations, the sigmoidal response of to becomes a hyperbolic curve, with a sharp, continuous increase of with at low . This change may reflect the necessity of remodelling chromatin – a process that requires ATP – for the transcription process. Chromatin remodelling has been noted by several studies [15], [16], [22] to play an important role in mRNA synthesis noise and hence downstream noise in gene expression. Rather than attempting to model this influence explicitly, we use the experimentally-determined form for to capture the overall dependence of transcription rate (including chromatin effects) on .

To summarise, in our model, transcription rate depends sigmoidally on ATP concentration – a relationship elucidated and quantified in recent experiments [31]. ATP concentration in turn depends linearly on the mitochondrial mass and functionality level of a cell and also on the cell volume. Cells with many, highly functional mitochondria will have higher levels of ATP and hence higher transcription rates than those with smaller, less functional mitochondrial populations.

Cellular dynamics

Our model for cell cycle dynamics consists of equations governing the time evolution of the key quantities volume, mitochondrial mass, and mitochondrial functionality. In the light of a recent study [53], and as cell cycle models often assume the exponential growth picture, we expect an exponential form for cell volume growth: . Here, is a function expressing the dependence of volume growth rate on other parameters.

We suggest that ATP concentration () plays a key role in powering growth of the cell, so cells with higher ATP levels have higher growth rates associated with cell volume and mitochondrial mass. This link postulates that biosynthesis rates are generally, like transcription, a function of ATP concentration. We note that although ATP concentration has been suggested [31] as a possible mechanism linking mitochondria and transcription rate, and some evidence supports this link, it may be the case that a different factor provides the causal mechanism, and ATP concentration is correlated with this underlying factor. For example, ROS, which adversely affect many cellular processes (including provoking a decrease in transcription rates [31]), may be an alternative to ATP, or a combination of ATP and ROS levels may act to determine transcription rate.

Numerous historical studies, both in HeLa cells [54] and other tissue types [55]–[60] have found that the density of mitochondrial mass (also called mitochondrial volume density) within cells of a given tissue type is consistent between generations and within populations. This consistency suggests that the time evolution of mitochondrial mass should be (a) coupled with the time evolution of volume and (b) of a form that allows damping of the inherent stochasticity at mitosis. In addition to these features, it is presumably reasonable to assume that mitochondrial growth is dependent on available (due to the required protein synthesis). We suggest a model that captures these required dependencies and incorporates mean-reversion, given by the dynamic equations:

(3)

(4)

where .

We note that this simple model does not distinguish between volume growth rates at different times in the cell cycle, but yields a smooth exponential growth in cell size throughout the cell cycle. We work in this picture for simplicity and generality, but note that a more sophisticated model would include a more detailed description of cell growth as another potential source of variability between cells.

The model's dynamics result (see Methods) in a convergence in mitochondrial density with time to a value .

Our simple model is sufficient to approximate a large set of experimental data

Here we list a set of comparisons between predictions from our model and experimental studies. Unless stated otherwise, we will use experimental data from the study of das Neves et al. [31], using the protocol ‘NX’ to refer to data in Fig. X of that study.

New experimental results are also consistent with this model

In Fig. 3, we also present new experimental results pertaining to our model. These new experiments were designed to characterise two additional features of cells in a population: a measure of the total level of mitochondrial function within cells and the modulation of cell cycle lengths by changing the oxidative state of the cell. The total level of mitochondrial function is experimentally measured using the intensity of signal from CMXRos, a dye that stains mitochondria and accumulates according to membrane-potential, integrated over a whole cell (see Methods). This signal reports on the integrated membrane potential across the entire cell, combining measures of mitochondrial mass and functionality. The population distribution of this quantity is of interest in exploring the link between mitochondrial mass and functionality between cells.

The modulation of cell cycle length with cellular oxidative state was investigated by observing the distribution of cell cycle lengths in a control population of cells and in populations of cells after anti-oxidant (dithiothretiol) or pro-oxidant (diamide) treatments (see Methods). Our model incorporates oxidative status by modulating the mean level of mitochondrial functionality, so mitochondria function more readily in an environment with low oxidative stress than one with high oxidative stress. As mitochondrial functionality is tied in our model, through growth rate, to cell cycle length, we would expect cell cycle lengths to decrease upon anti-oxidant treatment and increase upon pro-oxidant treatment.

Cell cycle lengths in different oxidative conditions

In Fig. 3J we show the mean and standard deviation of cell cycle lengths in a control population, and upon treatment with anti- and pro-oxidants (see Methods). In our simulations, these treatments are modelled by changing the value of , affecting the mean functionality of mitochondria (see Table 1). It is observed that treatment with anti-oxidants reduces cell cycle lengths, and treatment with pro-oxidants increases cell cycle lengths. In our model, this behaviour emerges from the dependence of the rate of volume growth on , and the increased levels resulting from mitochondria with higher functionality.

Table 1. Parameters and values employed in our model.

Parameter	Description	Value	Motivation
	memory term	0.5	Fit parameter – chosen to give a mean functionality of 1
	Sets mean functionality (control)	0.5	Fit parameter – chosen to give a mean functionality of 1
	Volume for mitosis (scale)		Fixed for consistency with maximum volume in Ref. [53]
	Proportionality between and		Fixed for consistency with mean ATP levels in Ref. [70]
	standard deviation at mitosis		Fixed by volume segregation data in Ref. [31]
	Set mean functionality (with anti-oxidant and pro-oxidant respectively)	(0.69, 0.09)	Fixed by transcription rate noise levels in Ref. [31]
	Fitting parameters for relationship between and	51.2, 44.7, ,	Fixed by functional form of in Ref. [31]
	growth rate		Chosen through optimisation – constrained by mean cell cycle length in Ref. [31].
	growth rate		Fixed ratio with through mitochondrial segregation data in Ref. [31]
	standard deviation at mitosis	0.34	Chosen through optimisation – constrained through transcription rate noise and cell cycle length variability in Ref. [31]

Open in a new tab

For further information see ‘Parameterisation of ’ and ‘Fitting Other Parameters’ in Text S1.

Summary of comparisons between experimental results and model predictions

It can be seen that several key experimental results require the inclusion of terms relating to mitochondrial variability for an explanation. In a situation without considerable mitochondrial influence on cellular variability, it may be expected that variability in cell cycle position among a population of unsynchronised cells may be a dominant source of noise. Physical distributions subject to such cell cycle noise would be expected to show a variance corresponding to an approximately twofold range, as this is the maximum difference in size between two unsynchronised cells. However, several results display data that varies over a considerably wider range than a factor of two, indicating that a factor other than cell cycle variability may be responsible. Most straightforwardly, Figs. 3A and 3I demonstrate pronounced cell-to-cell variability in the mass and functionality of mitochondrial populations. The distribution of transcription rate in Fig. 3I similarly shows a wide range of values.

Figs. 3D and 3E demonstrate the observed fact that mitochondrial inheritance at birth is a better predictor of cell cycle length than volume inheritance: an effect that relies on the presence of mitochondrial variability and mitochondrial influence on cellular growth. The variability in cell cycle length observed by modulating the oxidative state of the cell in Fig. 3J suggests that a source of variability that is sensitive to oxidative effects strongly affects cell cycle lengths. We believe that these results support the hypothesis that mitochondrial variability provides a significant contribution to the variability in distributions of the cellular properties we consider.

The correspondence between experimental data and the simulated behaviour of our model suggests that, although we have chosen simple functional forms in our model, the resulting behaviour is biologically relevant. However, we note here that our model was constructed from a phenomenological philosophy, with the intention of using experimental results to construct a plausible coarse-grained explanation for the influence of mitochondrial variability on extrinsic noise in general and transcription rate in particular: we were aware of all data from Ref. [31] when we were choosing the structure of our model though we only used a subset of available data to parameterise it. Our goal was to introduce a simplified but consistent mathematical summary of the data and to use this to motivate further experiments. To this end, we suggest a set of experiments in ‘Potential Experiments for Refinement’ (Text S1) that would support or contribute to further development of this model. We also note that many potential refinements could be made to our model and suggest several other functional forms in ‘Other Models’ (Text S1).

Noise in transcription rate depends on noise in mitochondrial segregation and functionality

We are now in a position to explore the dependence of the level of noise in transcription rate on the stochasticity in mitochondrial mass and function, and subsequent stochasticity in . To investigate the contribution of mitochondrial variability to transcription rate noise, we performed simulations of our model while varying , the variance associated with the inheritance of mitochondrial functionality, and , the variance associated with inheritance of mitochondrial mass. here gives the variance of the distribution by which mitochondrial mass is partitioned, and varying it under the assumption of binomial partitioning corresponds to changing the mitochondrial makeup of the cell: lower corresponds to more mitochondrial elements, each with smaller volume, while higher corresponds to fewer, larger mitochondrial elements, which are partitioned binomially at mitosis (see Methods).

In Fig. 5, the functional dependence of on mitochondrial variability ( and ) is shown from simulations. These results show that, for our model, the transcription rate noise is made up of significant contributions from both mitochondrial segregation and functionality. We also performed simulations where , the variability arising from uneven volume partitioning, was set to zero, and where cells were sampled at the same position in their cell cycle, removing different ages as a source of variability. As Fig. 5 shows, the removal of these sources of variability has little impact on the overall transcription rate noise level. These results lead us to suggest that mitochondrial sources of variability provide a strong contribution to cell-to-cell variability in transcription rate. This argument is supported by an approximate analytic treatment of the sources of error in transcription rate within our model (see ‘Estimating Noise Contributions’ in Text S1).

Mitochondrial variability can dominate noise in mRNA and protein expression

Having constructed and parameterised a model for mitochondrial variability and its effect on transcription in the cell, we now investigate the connection between these factors and downstream quantities: mRNA expression levels, and then (through further extension) protein expression levels. Noise in protein expression levels directly affects many cellular properties, as this noise causes cell-to-cell differences in the functional machinery available to perform cellular processes. Here we will investigate the influence of the mitochondrial variability suggested by the parameterisation of our model from experimental data on existing models for mRNA and protein expression. We connect our findings with the substantial existing body of literature on this topic in the Discussion section.

The production of mRNA and protein within a cell is often modelled using a master equation approach, addressing the probability of observing a given number of molecules at a given time. This analytical framework lends itself to the inclusion of our results for time-varying transcription rate (see Methods). Numerically, several studies have proposed techniques for incorporating time-varying rates in chemical kinetic systems [46], [47]: we use Shahrezaei et al.'s modification [45] to the Gillespie simulation method [63] to simulate our model system. This protocol allows us to investigate the relative importance of intrinsic contributions (resulting in differences in expression levels between identical genes within a single cell) and extrinsic contributions (resulting in differences in expression between identical genes in different cells in a population).

Fig. 6 shows the increase in mRNA expression (from a level of zero at the start of the simulation) from our analytic approach incorporating changing transcription rate, and in simulations run using (see Methods) a parameter set from Raj et al. [16], in two scenarios: one involving only intrinsic noise effects (no noise due to mitochondrial variability) and one involving extrinsic noise in transcription rate due to mitochondrial mass, functionality, and cell volume variability, of the magnitudes found through parameterising our model with experimental data. It can be seen that mitochondrial variability leads to a large increase in the total noise in mRNA expression levels: without extrinsic factors, the noise in mRNA expression at a given time ( hours) was , whereas with extrinsic factors. We note that the means for the intrinsic and extrinsic noise cases differ: this result is due to the nonlinear dependence of transcription rate on ATP concentration, so that distributions with the same mean but different variances may yield transcription rate distributions with different means.

We can also perform simulations on the more complicated system involving protein production (see ‘mRNA & Protein Levels’ in Text S1). With values from Raj et al. [16] for protein degradation and translation rate (see Methods), this approach allows us to simulate dual reporter experiments, where the expression of two distinct but identically regulated protein-encoding genes is measured. Each protein was translated from a different mRNA strand, so these simulations tracked four quantities: the expression levels of the two mRNAs and the two proteins. Simulations were performed on synchronised and asynchronous cells, and with set to their model values and set to zero. In these simulations, mRNA molecules and proteins were also distributed binomially between daughter cells at mitosis (see Methods).

Dual reporter simulations performed with the parameterisation chosen from Raj et al. [16] yield very low values for the magnitude of intrinsic noise. This low intrinsic noise was found to be due to the high copy number of proteins resulting from the parameterisation. To explore noise in systems with lower expression levels, we lowered the copy number of proteins by increasing the rates of mRNA and protein degradation (see Methods). Fig. 7 shows the resulting expression levels in two proteins with and without various sources of extrinsic noise, at the two different degradation rate protocols. These results show that, in our model, mitochondrial variability dominates the noise in protein expression levels. The spread of protein levels with mitochondrial and volume variability is much greater than the two-fold range achieved through cell cycle variability alone. Fig. 7 also illustrates that cells with higher mitochondrial mass and functionality generally have higher protein expression levels, though inheritance noise makes this correlation weaker.

In our model, we find that energy variability arising through mitochondrial stochasticity is the dominant source of variability in transcription rate, mRNA and protein expression levels. However, we note that the causal factors of stochasticity in mRNA and protein levels within the cell are significantly more complicated than the simple transcriptional model presented above. The rates of many of the processes involved in more extended models are functions of many factors which our model does not include. The inclusion of these complicating terms rapidly makes an analytic description of the model impossible. However, we note that stochastic simulation techniques may be used to explore the behaviour of complex model given estimates for the functional dependence of process rates on extrinsic variables [45].

We also note that several studies have observed a decrease in intrinsic noise at higher levels of protein expression [15], [27]. We do observe such a decrease, though in the default parameterisation the magnitude of this effect is very small owing to the consistently low intrinsic noise levels.

Mitochondrial noise, by modulating transcription rate, can affect stem cell differentiation

As an illustrative application of our model, demonstrating its physiological relevance, we consider how, through the extrinsic effects of [ATP] on protein levels, a link between mitochondrial content and stem cell differentiation behaviour may arise. Differentiation dynamics in stem cells have often been modelled as the result of expression asymmetries in lineage regulation genes that interact in a regulatory network [64]–[67], but the initial sources of this expression variability have not been clearly elucidated and are a topic of active debate. Here we show that transcription rate variability resulting from mitochondrial variability can affect the dynamics of expression of such control genes. Experimentally, a link between stem cell differentiation and mitochondria was suggested by a recent study in mouse embryonic stem cells [38], showing that pluripotent cells with low mitochondrial membrane potential had higher in vitro differentiation propensity, whereas those with higher membrane potential remained undifferentiated and formed large teratomas.

We explore two recent models for the cell fate decision between erythroid and myeloid cell fates directed by the cross-antagonistic master lineage regulators GATA1 and PU.1. One model, by Huang et al. [68], consists of a symmetric coupled ODE system for the expression levels of these two genes, including cross-repression and self-activation term (see Methods). Another model, by Chickarmane et al. [69], contains a similar but asymmetric ODE model, expanded to include interactions with a postulated third species which is promoted by GATA1 and represses PU.1. The Chickarmane et al. model also includes external signalling terms which may act to promote GATA1 and PU.1, and repress the third species. In these models, cell states are defined by the relative levels of expression of these genes, such that undifferentiated cells have comparable levels of each transcription factor, while the two differentiated cell types correspond to a state with high levels of one factor and low levels of the other. The interactions between the genes are parameterised by variables such as self-activation and cross-repression rates (see Methods). The phase space for both these models comprises three attractor basins, corresponding to the progenitor cell type and two differentiated cell types.

Within the Huang model, at low protein expression levels, smaller perturbations are required to shift attractor basins than at high expression levels – a feature consistent across a large range of parameterisations. Varying the parameterisation of the model (modelling differentiation-inducing signalling) changes the structure of these basins, so that the central undifferentiated basin becomes more or less stable to subsequent perturbation. We vary the default parameterisation of the model in an attempt to assess the effect of changes in transcription and translation rates (see Methods). We find that when the parameters related to the rate of production of proteins are low, the central, undifferentiated state is less stable than when they are high (see Fig. 8A), with a smaller volume of phase space leading to the undifferentiated basin.

Within the Chickarmane model, a different effect is observed. As before, we investigated the volume of phase space corresponding to the basin representing the undifferentiated state. We used a nonzero value for the external signalling term promoting PU.1 and explored the system at different transcription rates (see Methods). We found that increasing the transcription rate led to a decrease in the range of values of the external interaction which supported a stable undifferentiated state (see Fig. 8B). This decrease in the stability of the undifferentiated state arose from a smaller volume of phase space leading to the undifferentiated basin as transcription rate increased, with more phase space occupied by the GATA1 basin. This result contrasts with the increased stability of the Huang model at high transcription rate, due to the importance of the third species (the expression of which is dependent on transcription rate): at high transcription rate, the increased strength of the combined effect of self-activating GATA1 and production of the third species shifts the basin structure strongly towards GATA1.

These results suggest that cell-to-cell variability in mitochondrial mass and function may, through induced variability in transcription rate, have a significant effect on the stability of bipotent cells. If differentiation dynamics are asymmetric and involve an intermediate species (as in the Chickarmane model), we find that high transcription rate destabilises the undifferentiated state. This destabilisation may be viewed as a result of the increased sensitivity of the system to perturbations: the asymmetric regulatory architecture means that a small increase in GATA1 will be quickly amplified at high transcription rate, as more GATA1 and X are quickly produced. If differentiation dynamics are symmetric and do not involve another species (as in the Huang model), high transcription rates increase the width of the basin corresponding to the undifferentiated state, acting to stabilise this state. This stabilisation is due to the increased robustness to perturbations afforded by the high production rate of both species at high transcription rate: without asymmetric interactions, the higher expression level of both genes makes the system less responsive to small perturbations. The results that emerge from this symmetric case gives results that are qualitatively comparable to an experimental study [38] in which more cells with higher total mitochondrial membrane potential remained undifferentiated, suggesting that high mitochondrial performance stabilises the undifferentiated state.

Another, higher-order effect may conceivably play a role in both situations: several studies have found that, at high protein abundance levels (which may result from high transcription rates), intrinsic noise levels in protein expression decrease. While the parameterisation of our dual reporter studies is such that these effects are small, the fact that less noise is expected at higher protein expression levels suggests a third mechanism by which high mitochondrial content may stabilise pluripotent cells. The contrasting results highlight the potential of experimental investigation of the effects of global transcription rate on the stability of multipotent states to inform of additional qualitative behaviors that models of lineage decision should be expected to exhibit.

Conclusions

We find, through a phenomenological model constructed to reproduce recent data on mitochondrial and ATP variability, that stochastic inheritance of mitochondria at mitosis and variability in mitochondrial function may be important sources of noise in transcription. By extension, these factors may contribute significantly to noise in protein expression further downstream. We have proposed experimental tests to refine our model and demonstrate its application in existing models for mRNA and protein production and stem cell differentation, and discussed how these findings integrate into the current understanding of extrinsic noise in cellular biology. In particular, what our paper suggests is the need for multimodal single cell experiments through time (and through division) investigating coarse-grained measures of energy status, cellular volume, mitochondrial mass, and global rates of transcription and translation. Cellular variability is of central physiological importance but we suggest that to understand this we must elucidate the relationships between certain core variables, including the relationship between the machinery of expression and degradation and the energy status of the cell.

Parameterisation of our model from experimental data

We used a subset of available experimental data (shown in Fig. 2) to choose numerical values for our key parameters. Some parameter values were fixed, in the sense that maximising the agreement between predictions from our model and a single experimental study allowed an optimal value to be chosen straightforwardly ( and the ratio were parameterised by data from Ref. [31] (respectively the sigmoidal relationship between and transcription rate, the variance associated with volume partitioning, and the variance associated with mitochondrial partitioning), by the maximal cell volume observed in data from Ref. [53], and by comparing the mean properties of simulated cells to mean ATP levels reported in a population of HeLa cells in Ref. [70]). Other parameters influenced a range of predictions and values for these parameters were chosen by optimising the fit to experimental data across the set of results that they influenced (see ‘Fitting Other Parameters’ in Text S1).

Table 1 summarises the parameters and values employed in our model.

Model specifics

In the following we provide details of our model. The key equations specifying the key dynamics of cells in our model are given by Eqns. 1–4. The coupled ODEs in our model admit an analytic solution by simple integration:

(5)

(6)

The reader will note that a variety of models for mitochondrial and volume growth will yield similar forms. In addition, other functional forms for the level of ATP may be suggested. A selection of these alternative forms are explored in ‘Other Models’ in Text S1. Within this model, at mitosis,

(7)

(8)

where is a normal distribution with mean and variance , and determine the daughter volumes and mitochondrial masses . The variances of the mitochondrial distribution is chosen to represent a binomial distribution (, with , for high ). The variance of the volume distribution is chosen to match experimental data on volume partitioning.

Functionality evolves through changes at mitosis events, which follow an AR(1) process. A daughter cell's functionality is determined from its parent's functionality :

(9)

Upon mitosis, both daughter cells inherit the same value, drawn from Eqn. 9. Once chosen, a cell's functionality remains constant throughout one cell cycle. We can model treatment with anti- and pro-oxidants by respectively increasing and decreasing , raising or lowering the mean functionality of mitochondria. With , there is a very small but finite chance that cells will inherit zero (or lower) functionality. To avoid this unphysical case, we impose a cutoff of on mitochondrial functionality. Values under this cutoff are resampled from Eqn. 9. Another model, where varies continuously within a cell cycle, is considered in ‘Other Models’ in Text S1.

Virtual mitochondria

In our model, the standard deviation in mitochondrial mass at mitosis is a function of the number of virtual mitochondria in the cell, , from binomial partitioning. We can vary the number of virtual mitochondria in the cell while keeping the total mitochondrial volume constant, by changing the volume assigned to an individual virtual mitochondrion. Let an individual virtual mitochondria have volume , and let there be virtual mitochondria in a cell. The total mitochondrial volume is , and the mean inherited mitochondrial volume is half of this. The standard deviation in virtual mitochondrial number from a binomial distribution of the virtual mitochondria is , so the standard deviation of mitochondrial volume is . A standard deviation of then corresponds to virtual mitochondria per cell.

Cell dynamics simulations

To simulate a population of cells, we used a simple Euler method to solve the dynamic equations. A population of cells was simulated. When mitosis occurred, a random cell from the population was chosen for replacement by the new cell. To measure distributions, this procedure was continued until the distributions stabilised. To measure sister-sister and between-generation correlations, a list of relationships between cells was maintained, with sister pairs only sampled when both sisters underwent an entire cell cycle without replacement. We note that this removal of randomly-chosen cells may potentially introduce artefacts into the results, as the constant probability of cell removal means that the probability of a cell surviving to a certain age is decreasing. We checked our results with a different simulation protocol: running the system and allowing exponential growth up to cells, with no removal. The resulting distributions were indistinguishable from the first protocol, showing that the removal rate is low enough so that the statistics are not affected.

mRNA & protein expression simulations

For the simple system illustrating mRNA expression levels alone, cells were simulated using the modification of Shahrezaei et al.'s [45] to the Gillespie simulation method [63] in two scenarios: one involving only intrinsic noise () and one involving extrinsic noise due to mitochondrial mass, functionality, and cell volume variability. In both cases, the initial copy number of mRNA molecules was set to zero, to illustrate differences in the dynamics of the system: transcription started at the birth of the cell, and there were no effects from mRNA inheritance. In the intrinsic noise experiments, a population of cells were simulated from identical initial conditions, with their values set to the means of those variables obtained from simulations. The variability in mRNA expression levels in these cells was therefore solely due to intrinsic noise. For the extrinsic noise simulations, mRNA expression was simulated in the heterogeneous population of cells that resulted from dynamic simulation. An ensemble of cells was analysed for both cases and the mRNA content at each timestep recorded.

To simulate the more complicated systems and investigate inheritance effects, we coupled the Shahrezaei et al. simulation protocol with the simple ODE solver so that simulation of a population of cells growing and producing mRNA and protein involved the following algorithm. 1) Use the ODE solver to calculate a cell's time of mitosis and the time series of volume and transcription rate throughout the cell lifetime. 2) Use Shahrezaei et al.'s method to compute the time behaviour of mRNA and protein levels given these time series for production rates. 3) Create daughter cells with noisy partitioning of volume, mitochondrial mass (binomial) and functionality (AR(1)), and mRNA and protein copy numbers (binomial).

After Raj et al. [16], we employ the following model values for birth () and death () rates of mRNA () and protein (): . In the case of birth rates, we used these mean values to scale the curves from das Neves et al. [31] so that the mean level observed in a population gave the mean values from Raj et al. (see ‘mRNA & Protein Levels’ in Text S1). We also ran experiments with the degradation rates increased 100-fold: , to explore the behaviour of the system at lower expression levels.

Experiments

CMXRos labelling was done according to manufacturer (invitrogen) instructions, cells were incubated with the probe (75 nM) for 15 min. at 37 C and washed with warm PBS.

For the dual dye experiments, cells were loaded simultaneously with CMXRos and MitoTracker Green FM dye (Molecular Probes) for 20 min and after a brief PBS rinse fixed in PBS with 4% paraformaldehyde.

Cell cycle length experiments were done by growing cells in the presence of media alone or containing either the anti-oxidant Dithiothreitol (250 microM DTT) or the pro-oxidant Diamide (50 microM). Cell cycle length was measured as the time interval between two mitotic events in a single cell, analysed by live cell imaging using the Cell IQ platform and image analysis software (Chipman Tech.).

In the flow cytometry experiments, tripsinised Hela cells were washed and incubated with 20 nM MitoTracker Green FM dye (Molecular Probes) in medium at 37 C for 15 min, and then washed. Cells were analyzed using flow cytometry (Dako CyAn ADP).

Master equations

Let us consider the master equation for the transcription process, describing the probability of observing the system with mRNAs at time :

(10)

where is the probability of observing mRNAs at a given time , is transcription rate, and is an mRNA degradation rate.

Using a linear approximation for (see ‘Parameterisation of ’ in Text S1), we can solve this by using a generating function approach (see ‘mRNA & Protein Levels’ in Text S1). The mean, variance, and probability distribution of mRNA copy number at arbitrary time are given by:

(11)

(12)

(13)

where is the Kummer confluent hypergeometric function, and .

Stem cell model

The progenitor cell differentation model of Huang et al. [68] consists of the following equations for the evolution of protein expression levels (GATA1) and (PU.1):

(14)

(15)

In this parameter set, variables are self-activation rates, variables are cross-repression rates, are decay rates, and and control the functional form of these processes. Huang et al. show that altering these parameters changes the structure of the corresponding attractor landscape, so that the central undifferentiated attractor basin changes in size, affecting the predisposition of the system to differentiate. This landscape change gives a ‘priming’ of the system such that the effect of subsequent asymmetries may vary. In this study, we vary the landscape by symmetrically varying and , the parameters associated with activation and repression, from the default parameterisation for . Specifically, we modulated and with a multiplicative factor : . We draw the connection between higher values of and and higher transcription and translation rates, as the rate of production of chemical species is increased by an increase in these parameters.

The model of Chickarmane et al. [69] involves a relationshup between three dynamic variables:

(16)

(17)

(18)

where are the concentrations of GATA1, PU.1 and a postulated chemical species X respectively. are external signalling factors: and promote GATA1 and PU.1 respectively, and represses X. Other variables take default values , , , , , . To vary transcription rates, we modulate and terms with a multiplicative factor which we take to be proportional to transcription rate . We used with , allowing an external signal that promotes PU.1, to promote stability of the undifferentiated state.