Stochastic modelling, Bayesian inference, and new in vivo measurements elucidate the debated mtDNA bottleneck mechanism

Dangerous damage to mitochondrial DNA (mtDNA) can be ameliorated during mammalian development through a highly debated mechanism called the mtDNA bottleneck. Uncertainty surrounding this process limits our ability to address inherited mtDNA diseases. We produce a new, physically motivated, generalisable theoretical model for mtDNA populations during development, allowing the first statistical comparison of proposed bottleneck mechanisms. Using approximate Bayesian computation and mouse data, we find most statistical support for a combination of binomial partitioning of mtDNAs at cell divisions and random mtDNA turnover, meaning that the debated exact magnitude of mtDNA copy number depletion is flexible. New experimental measurements from a wild-derived mtDNA pairing in mice confirm the theoretical predictions of this model. We analytically solve a mathematical description of this mechanism, computing probabilities of mtDNA disease onset, efficacy of clinical sampling strategies, and effects of potential dynamic interventions, thus developing a quantitative and experimentally-supported stochastic theory of the bottleneck. DOI: http://dx.doi.org/10.7554/eLife.07464.001

Dangerous damage to mitochondrial DNA (mtDNA) can be ameliorated during mammalian development through a highly debated mechanism called the mtDNA bottleneck.Uncertainty surrounding this process limits our ability to address inherited mtDNA diseases.We produce a new, physically motivated, generalisable theoretical model for mtDNA populations during development, allowing the first statistical comparison of proposed bottleneck mechanisms.Using approximate Bayesian computation and mouse data, we find most statistical support for a combination of binomial partitioning of mtDNAs at cell divisions and random mtDNA turnover, meaning that the debated exact magnitude of mtDNA copy number depletion is flexible.New experimental measurements from a wild-derived mtDNA pairing in mice confirm the theoretical predictions of this model.We analytically solve a mathematical description of this mechanism, computing probabilities of mtDNA disease onset, efficacy of clinical sampling strategies, and effects of potential dynamic interventions, thus developing a quantitative and experimentallysupported stochastic theory of the bottleneck.
Mitochondria are vital energy-producing organelles within eukaryotic cells, possessing genomes (mitochondrial DNA or mtDNA) that replicate, degrade and develop mutations [1,2].MtDNA mutations have been implicated in numerous pathologies including fatal inherited diseases and ageing [3,4,5,1].Combatting the buildup of mtDNA mutations is of paramount importance in ensuring an organism's survival.Substantial recent medical, experimental and media attention focused on methods to remove [6] or prevent the inheritance of [7,8,5,9,10] mutated mtDNA in humans.
One means by which organisms may ameliorate the mtDNA damage that builds up through a lifetime is through a developmental process known as bottlenecking.Immediately after fertilisation, a single oocyte (which may contain > 10 5 individual mtDNAs) may have a nonzero mtDNA mutant load or heteroplasmy (the proportion of mutant mtDNA in the cell).As the number of cells in the developing organism in-creases, the intercellular population will then acquire an associated heteroplasmy variance, that is, the variance in mutant load across the population of cells (Fig. 1A), allowing removal of cells with high heteroplasmy and retention of cells with low heteroplasmy.Intense and sustained debate exists as to the mechanism by which this increase of heteroplasmy variance occurs.Several experimental results in mice suggest that, during development, the copy number of mtDNA per cell in the germ cell line drops dramatically to ∼ 10 2 , reducing the effective population size of mitochondrial genomes [11,12].One postulated bottlenecking mechanism is that this low population size accelerates genetic drift and so increases the cell-to-cell heteroplasmy variance [11,13,14,15], which was first observed to generally increase from primordial germ cells through primary oocytes to mature oocytes [16].However, independent experimental evidence [12] suggests that heteroplasmy variance increases negligibly during this copy number reduction, though this interpretation has been debated [17].Ref. [12] shows heteroplasmy variance rising during folliculogenesis, after the mtDNA copy number minimum has been passed.In yet another picture, supported by conflicting experimental results [18,19], heteroplasmy variance increases with a less pronounced decrease in mtDNA copy number (a minimum copy number > 10 3 in mice), solely through random effects associated with partitioning at cell divisions.Clearly a consensus on this important mechanism is yet to be reached.
Important existing theoretical work on modelling the bottleneck has assumed a particular underlying mechanism [13,20] or derived statistics of mtDNA populations [21,14,22,23] without explicitly considering changing mtDNA population size, or the discrete nature of the mtDNA population: effects which may powerfully affect mtDNA statistics.To capture these effects it is necessary to employ a 'bottom-up' physical description of mtDNA as populations of individual, discrete elements subject to replication and degradation, as in, for example, Refs.[21] and [24].Exploring the bottleneck also requires explicitly modelling partitioning dynamics throughout 1 arXiv:1512.02988v1[q-bio.QM] 9 Dec 2015 a series of cell divisions, over which population size can change dramatically.While previous simulation work [11,25] has taken such a philosophy with specific model assumptions, we are not aware of such a study allowing for the wide variety of replication and partitioning dynamics proposed in the literature; we further not that replication-degradation-partitioning mtDNA models are yet to be fully described analytically.Nor is there a general quantitative framework under which different proposed bottleneck mechanisms can be statistically compared given extant data (although statistical analyses focusing on particular mechanisms and individual sets of experimental results have been used throughout the literature, for example, using a Bayesian approach under a particular bottleneck model to infer model bottleneck size [26]).Combined developments in theory and inference are therefore required to make progress on this important question.
We remedy this situation by constructing a general model (features and parameters described in Fig. 1) for the population dynamics of the bottleneck, able to describe the range of proposed mechanisms existing in the literature.Using experimental data on mtDNA statistics through development [16,11,12,18], we use approximate Bayesian computation [27,28,29,30] to rigorously explore the statistical support for each mechanism, showing that random mtDNA turnover coupled with binomial partitioning of mtDNAs at cell divisions is highly likely, and that the debated magnitude of mtDNA copy number reduction is somewhat flexible.Subsequently, we confirm the predictions of this model by performing new experimental measurements of heteroplasmy statistics in mice with an mtDNA admixture, including a wild-derived haplotype, that is genetically distinct from previous studies.We then analytically solve the equations describing mtDNA population dynamics under this mechanism and show that these results allow us to investigate potential interventions to modulate the bottleneck (suggesting that upregulation of mtDNA degradation may increase the power of the bottleneck to avoid inherited disease; we discuss potential strategies for such an intervention) and yield quantitative results for clinical questions including the timescales and probabilities of disease onset, and the efficacy of strategies to sample heteroplasmy in clinical planning.

A general mathematical model encompassing proposed bottlenecking mechanisms
We will consider three different classes of proposed generating mechanisms for the mtDNA mechanisms: namely, those proposed in Cree et al. [11]; Cao et al. [18] and Wai et al. [12].We will refer to these mechanisms by their leading author name.The Cree mechanism involves random replication and degradation of mtDNAs throughout development, and binomial partitioning of mtDNAs at cell divisions.The Cao mechanism involves partitioning of clusters of mtDNA at each cell division, thus providing strong stochastic effects associated with each division.We consider a general set of dynamics through which this cluster inheritance may be manifest, including the possibility of heteroplasmic 'nucleoids' of constant internal structure [31], sets of molecules or nucleoids within an organelle, homoplasmic clusters, and different possible cluster sizes (see Supplementary Information).The Wai mechanism involves the replication of a subset of mtDNAs during folliculogenesis.We note that this latter mechanism can be manifest in several ways: (a) through slow random replication of mtDNAs (so that, in any given time window, only a subset of mtDNAs will be actively replicating) or (b) through the restriction of replication to a specific subset of mtDNAs at some point during development.We will refer to these different manifestations as Wai (a) and Wai (b) respectively.The Wai (a) mechanism and the Cree model can both be addressed in the same mechanistic framework (with potentially different parameterisations): if the rate of random replication in the Cree model is sufficiently low during folliculogenesis, only a subset of mtDNAs will be actively replicating at any given time during this period, thus recapitulating the Wai (a) mechanism (see Supplementary Information).We will henceforth combine discussion of the Wai (a) and Cree mechanisms into what we term the birth-death-partition (BDP) mechanism.
We seek a physically motivated mathematical model for the bottleneck that is capable of reproducing each of these mechanisms.Our general model for the bottleneck (detailed description in Methods) involves a 'bottom-up' representation of mtDNAs as individual intracellular elements capable of replication and degradation (Fig. 1B) with rates λ and ν respectively.A parameter S determines whether these processes are deterministic (specific rates of proliferation) or stochastic (replication and degradation of each mtDNA is a random event).These rates of replication and degradation of mtDNA are likely strongly linked to mitochondrial dynamics within cells, through the action of mitochondrial quality control [32,33] modulated by mitochondrial fission and fusion [34,35,36], which can act to recycle weaklyperfoming mitochondria [37,38].This quality control can be represented through the degradation rates assigned to each mtDNA species, which may differ (for selective quality control) or be identical (for non-selective turnover).
The proportion of mtDNAs capable of replication is controlled by a parameter α in our model, dictating the proportion of mtDNAs that may replicate after a cutoff time T .Thus, if α = 1, all mtDNAs may replicate; if α < 1, replication of a subset proportion α of mtDNAs is enforced at this cutoff time.At cell divisions, mtDNAs may be partitioned either deterministically, binomially, or in clusters according to a parameter c (Fig. 1C).
The copy number of mtDNA per cell is observed to vary dramatically during development, with dynamic phases of copy number depletion and different rates of subsequent recovery observed.Additionally, cell divisions occur in the germline at different rates during development, with cells becoming largely quiescent after primary oocytes develop.To explicitly model these different dynamic regimes, and the behaviour of mtDNA copy number during each, we include six different dynamic phases throughout development, MtDNA copy number (j = 1 wildtype; j = 2 mutant) MtDNA replication; degradation rate in phase i Length of a cell cycle in phase i (if i is a cycling phase) Describes mtDNA dynamics: 0 (deterministic replication and degradation); 1 (stochastic replication and degradation) Number of cell cycles in phase i (0 if i is a quiescent phase) Describes mtDNA partitioning at cell divisions: 0 (deterministic); >0 (binomial partitioning of clusters of size c) Initial heteroplasmy; initial mtDNA copy number (at conception) The proportion of mtDNAs capable of replication after time T m j h P(x); E(x); V(x); V'(x)  each with different rates of replication and degradation (labelled with subscript i labelling the dynamic phase: hence λ 1 , ν 1 , ..., λ 6 , ν 6 ), and allowing for different rates of cell division or quiescence.This protocol enables us to explicitly model effects of changing population size throughout development rather than assuming dependence on a single, coarsegrained effective population size; and to include the effects of specific and varying cell doubling times.A summary of symbols used in our model and throughout this article is presented in Fig. 1D.
Our model, with suitable parameterisation, can thus mirror the dynamics of the Cree and Wai (a) mechanisms (stochastic dynamics and binomial partitioning, which we refer to as the BDP mechanism); the Cao mechanism (clustered partitioning); and Wai (b) mechanism (deterministic dynamics, restricted subset of replicating mtDNAs).The Cao mechanism, partitioning of clusters of mtDNA molecules, represents the expected case if mtDNA is partitioned in colocalised 'nucleoids' within each organelle (or in other sub-organellar groupings).The size of mtDNA nucleoids is debated in the literature [1,39,40] (although recent evidence from highresolution microscopy suggests that nucleoid size is generally < 2 [41], consonant with recent evidence that individual nucleoids may be homoplasmic [42]); our model allows for inheritance of homoplasmic or heteroplasmic nucleoids of arbitrary characteristic size c, thus allowing for a range of suborganellar mtDNA structure.We discuss the impact of mixed or fixed nucleoid content in the Supplementary Information.

A birth-death-partition model of mtDNA dynamics has most statistical support given experimental measurements
We take data on mtDNA copy number in germ line cells in mice from three recent experimental studies [11,18,12].We also use data from two experimental studies on heteroplasmy variance in the mouse germ line during development [16,12].This data, by convention [17], is normalised by heteroplasmy level h, giving where normalised variance V (h) is a quantity that will be often used subsequently.This normalised variance controls for the effect of different or changing mean heteroplasmy, and thus allows a comparison of heteroplasmy variance among samples with different mean heteroplasmies and subject to heteroplasmy change with time.We use a time of 100 dpc to correspond to mature oocytes (see Methods).These heteroplasmy variance studies employ intracellular combinations of the same pairing of mtDNA haplotypes (NZB and BALB/c), modelling two different mtDNA types within a cellular population.We take data on cell doubling times from a classical study [43] (see Methods).A possible summary of these data (although they provoke ongoing debate; see Discussion) is that, as shown in Fig. 2A, the existing data on normalised heteroplasmy variance shows initially low variance until ∼ 7.5 dpc (days post conception, which we use as a unit of time throughout), rising to intermediate values between 7.5 and 21 dpc, gradually rising further subsequently to become large  in the mature oocytes of the next generation.In Fig. 2A, and throughout this article, experimentally measured data will be depicted as circular points and inferred theoretical behaviour will be depicted as lines or shaded regions.
Fig. 2A shows mtDNA population dynamic trajectories resulting from optimised parameterisations of each of the mechanisms we consider (see Methods).In Fig 2B we show posterior probabilities on each of these mechanisms.These posterior probabilities give the inferred statistical support for each mechanism, derived from model selection performed with approximate Bayesian computation (ABC) [27,28,29,30] using uniform priors.ABC involves choosing a threshold value dictating how close a fit to experimental data is required to accept a particular model parameterisation as reasonable.In our case, this goodness-of-fit is computed using a comparison of squared residuals associated with the trajectories of mean mtDNA copy number and normalised heteroplasmy variance (see Methods and Supplementary Information).Each of the experimental measurements corresponds to a sample variance, derived from a finite number of samples of an underlying distribution of heteroplasmies, and therefore has an associated uncertainty and sampling error [14].The reasonably small sample sizes used in these sample variance measurements are likely to underestimate the underlying heteroplasmy variance (the target of our inference).Our ABC approach naturally addresses these uncertainties by using summary statistics derived from sampling a set of stochastic incarnations of a given model, where the size of this set is equal to the number of measurements contributing to the experimentally-determined statistic (see Methods).Fig. 2B clearly shows that as the ABC threshold is decreased, requiring closer agreement between the distributions of simulated and experimental data, the posterior probability of the BDP model increases, to dramatically exceed those of the other models.This increase indicates that the BDP model is the most statistically supported, and capable of providing the best explanation of experimental data (which can be inutitively seen from the trajectories in Fig. 2A).We note that ABC model selection automatically takes model complexity into account, and conclude that the BDP mechanism is the best supported proposed mechanism for the bottleneck.Briefly, this result arises because the BDP model produces increasing variance both due to early cell division stochasticity and later random turnover.By contrast, the Cao model alone only increases variance in early development when cell divisions are occurring.Qualitatively, this behaviour through time holds regardless of cluster (nucleoid) size and regardless of whether clusters are heteroplasmic or homoplasmic (allowing heteroplasmic clusters decreases the magnitude of heteroplasmy variance but not its behaviour through time, see Supplementary Information).The Wai (b) model alone similarly only increases variance at a single time point (later, during folliculogenesis).
In Ref. [12], visualisations of cells after BrU incorporation show that a subset of mitochondria retain BrU labelling, which the authors suggest indicates that a subset of mtD-NAs are replicating.In the Supplementary Information, we show that the BDP model also results in the observation of only a subset of replicating mtDNAs over the timeframe corresponding to these experimental results.These observations thus correspond to results expected from the random turnover from the BDP model.We also note the mathematical observation that the Wai (b) mechanism requires the replication of < 1% of mtDNAs during folliculogenesis to yield reasonable heteroplasmy variance increases (Fig. 2A shows the optimal case with α = 0.006; optimal fits to data generally show 0.005 < α < 0.01), and the proportions of loci visible in Ref. [12] are substantially higher than this required 1% value.
We show in the Supplementary Information that the heteroplasmy statistics corresponding to binomial partitioning also describe the case where the elements of inheritance are heteroplasmic clusters, where the mtDNA content of each cluster is randomly sampled from the population of the cell (either once, as an initial step, or repeatedly at each division).This similarity holds broadly, regardless of whether the internal structure of clusters is constant across cell divisions or allowed to mix between divisions.The BDP model, in addition to describing the partitioning of individual mtDNAs, also thus represents the statistics of mtDNA populations in which heteroplasmic nucleoids are inherited [31], or individual organelles containing a mixed set of mtDNAs or nucleoids are inherited, regardless of the size of these nucleoids (see Discussion).

Parameterisation and interpretation of the birth-death-partition model
Having used ABC model selection to identify the BDP model as the most statistically supported, we can also use ABC to infer the values of the governing parameters of this model given experimental data.Figs.A and B shows the trajectories of mean copy number and mean heteroplasmy variance resulting from model parameterisations identified through this process.Fig. C shows the inferred behaviour of mtDNA degradation rate ν in the model, a proxy for mtDNA turnover (as the copy number is constrained).Turnover is generally low during cell divisions, allowing heteroplasmy variance to increase due to stochastic partitioning.Turnover then increases later in germ line development, resulting in a gradual increase of heteroplasmy variance after birth until the mature oocytes form in the next generation.
Fig. D shows posterior distributions on the copy number minimum and total turnover (see Methods) resulting from this process; posteriors on all other parameters are shown in Supplementary Information.Substantial flexibility exists in the magnitude of the copy number minimum, illustrating that observed heteroplasmy variance can result from a range of bottleneck sizes from ∼ 200 to > 10 3 ; going some way towards reconciling the conflict between Refs.[18] and [19] and Refs.[11] and [12].The total amount of mtDNA turnover (presented as σ = 6 i=3 ν i τ i , the product of turnover rate and the time for which this rate applies, summed over quiescent dynamic phases; for example, a turnover rate of 0.1 hr −1 for 30 days yields σ = 0.1 × 24 × 30 = 72) is constrained more  than the specific trajectory of mtDNA turnover rates, showing that a variety of time behaviours of turnover are capable of producing the observed heteroplasmy behaviour.

Experimental verification of the birth-deathpartition model
The bottleneck mechanism identified through our analysis has several characteristic features which facilitate experimental verification.Key among these are the prediction that heteroplasmy variance acquires an intermediate (nonzero, but not maximal) value as a result of the copy number bottleneck, then continues to increase due to mtDNA turnover in later development.Our theory also produces quantitative predictions regarding the structure of heteroplasmy distributions at arbitrary times.
The existing data that we used to perform inference and model selection display a degree of internal heterogeneity, coming from several different experimental groups.Furthermore, these data represent statistics resulting from a single pairing of mtDNA types, and it is thus arguable how conclusions drawn from them may represent the more genetically diverse reality of biology.Ref. [44] recently addressed the issue of this limited number of mtDNA pairings by producing novel mouse models involving mixtures of standard and several new, unexplored, wild-derived haplotypes which capture a range of genetic diversity.To test the applicability and generality of our predictions, we have perfomed new experimental measurements of germline heteroplasmy variance in these model animals under a consistent experimental protocol (see Methods).We use the 'HB' mouse line from Ref. [44] pairing a wild-derived mtDNA haplotype (labelled 'HB' after its source in Hohenberg, Germany) with C57BL/6N; we refer to this model as 'HB'.
Heteroplasmy measurements were taken in oocytes sampled from mice at ages 24-61dpc (see Methods and Supplementary Information; raw data in Source data 1).The statistics of these measurements yielded E(h), V(h) and V (h) as previously.This age range was chosen to address the regions with most power to discriminate between the competing models; the existing V (h) data is most heterogeneous around 20-30dpc and the later datapoints allow us to detect developmental heteroplasmy behaviour after the copy number minimum.Fig.A shows these V (h) measurements.The qualitative behaviour predicted by the BDP mechanism is clearly visible: variance around birth (after the copy number bottleneck) is low but non-zero, subsequently increasing with time.The ability of the BDP model to account for the magnitudes and time behaviour of heteroplasmy variance more satisfactorily than the alternative models is shown by the model fits in Fig.
A. We explored these new data quantitatively through the same model selection approach used for the existing data.As shown in Fig. B, the BDP mechanism again experiences by far the strongest statistical support in this genetically different system.Fig. C shows the result of our parameteric inference approach using these V (h) measurements coupled with the E(m) measurements used previously (employing our assumption that modulation of copy number by heteroplasmy in this non-pathological haplotype is small).Strikingly, the quantitative behaviour of V (h) with time inferred from the HB model (red) matches the previous behaviour inferred from the NZB/BALB/c system (blue) very well, suggesting that our theory is applicable across a range of genetically distinct pairings.We note that the shaded region in Fig. C corresponds to credibility intervals around the mean behaviour of V (h), and the fact that individual V (h) datapoints (subject to fluctuations and sampling effects) do not all lie within these intervals is not a signal of poor model choice.An analogous situation is the observation of a scatter of datapoints outside the range of the standard error on the mean (s.e.m.), which does not imply a mistake in the s.e.m. estimate.The difference between the trace in shows the distribution of model behaviours over the posterior distributions on parameters: the mean V (h) trace of this distribution is comparable but not equivalent to that from the single best-fit parameterisation.
To confirm more detailed predictions of our model, we also examined the specific distributions of heteroplasmy in our new measurements.Given a mean heteroplasmy and an organismal age, the parameterised BDP model predicts the structure of the heteroplasmy distribution (see Methods and  [44] is used to compare distributions with different mean heteroplasmy.Red jitter points are samples from sets used to parameterise the BDP model; red curvesshow the 95 % range on transformed heteroplasmy with time inferred from these samples.Blue jitter points are samples withheld independent from this parameterisation; their distributuions fall within the independently inferred range.Insets show, in untransformed space, distributions of the withheld heteroplasmy measurements (blue) compared to parameterised predictions (red); no withheld datasets show significant support against the predicted distribution (Anderson-Darling test, p < 0.05).
next section).We parameterised the model using V (h) values from a subset of half of the new measurements (chosen by omitting every other sampled set when ordered by time).Fig. D shows a comparison of measured heteroplasmy distributions with a 95 % bound from the parameterised BDP model.We then tested the predictions of the parameterised model against the other half of new measurements.8 of the test measurements (2.4%) fell outside the inferred 95% bound from the training dataset, illustrating a good agreement with distributional predictions.The Anderson-Darling test was used to compare the distribution of heteroplasmy in sampled oocytes with distributions predicted by our theory (given age and mean heteroplasmy); no set of samples showed significant (p < 0.05) departures from the hypothesis that the two distributions were identical.Some example distributions are presented in Fig.

The birth-death-partition model is analytically tractable
Importantly, the birth-death-partitioning model yields analytic solutions for the values of all genetic properties of interest, using tools from stochastic processes (detail in Methods and Supplementary Information).These results facilitate straightforward further study and fast predictions of timescales and probabilities of interest.The full theoretical approach is detailed in the Supplementary Information, and equations for the mean and variance of mtDNA populations and heteroplasmy are given in the Methods.In Fig. 2A we illustrate that these analytic results exactly match the numeric results of stochastic simulation, a result that holds across all BDP model parameterisations.It is also straightforward to calculate the fixation probability P(m = 0), which allows us to characterise all heteroplasmy distributions that arise from the bottlenecking process, even when highly skewed (see Methods and Supplementary Information).We have thus obtained analytic solutions for the time behaviour of mtDNA copy number and heteroplasmy throughout the bottleneck with no assumptions of continuous population densities or fixed population size, under a physical model with the most statistical support given experimental data.

Mitochondrial turnover, degradation, and selective pressures exert quantifiable influence on heteroplasmy variance
We can use our theory to explore the dependence of bottleneck dynamics on specific biological parameters.We first explore the effects of modulating mtDNA turnover by varying λ and ν in concert, corresponding to an increase in mtDNA degradation balanced by a corresponding increase in mtDNA replication.This increased mtDNA turnover increases the heteroplasmy variance due to bottlenecking (see Fig. 4A).This result arises due to the increased variability in mtDNA copy number from the underlying random processes occurring at increased rates.Additionally, we find that increasing mtDNA degradation ν without increasing λ also increases heteroplasmy variance, in addition to decreasing the overall mtDNA copy number (Fig. 4B).Applying this unbalanced increase in mtDNA degradation without a matching change in replication has a strong effect on mtDNA dynamics as it corresponds to a universal change in the 'control' applied to the system, analogous, for example, to changing target copy numbers in manifestations of relaxed replication [21].The simple model we use does not include feedback and controls mtDNA dynamics solely through kinetic parameters.Perturbing the balance of these parameters thus strongly affects the expected behaviour of the system.As we discuss later, elucidation of the specific mechanisms by which control is manifest in mtDNA populations will require further research, but these numerical experiments attempt to represent the cases where a perturbation is naturally compensated for (matched changes, Fig. 4A) and where it is not (unbalanced change, Fig. 4B).
These results suggest that an artificial intervention increasing mitochondrial degradation may generally be expected to increase heteroplasmy variance during development.An increase in mtDNA degradation is expected to either directly increase heteroplasmy variance (Fig. 4B) if mtDNA populations are weakly controlled, or to provoke a compensatory, population-maintaining increase in mtDNA replication, thus increasing mtDNA turnover, which also acts to increase variance (Fig. 4C) if mtDNA populations are subject to feedback control.The increase in variance through either of these pathways will increase the power of cell-level selection to remove cells with high heteroplasmy and thus purify the population.For this reason, we speculate that mitochondrial degradation may represent a potential clinical target to address the inheritance of mtDNA disease (more detail in Supplementary Information).
Our model also allows us to explore the effect of different mtDNA types experiencing different selective pressures, by setting λ 1 = λ 2 (mutant mtDNA experiences a proliferative advantage or disadvantage).Such a selective difference causes changes in both mean heteroplasmy and heteroplasmy variance, as shown in Fig. 4C (for example, if heteroplasmy decreases towards zero, heteroplasmy variance will also decrease, as the wild type is increasingly likely to become fixed).We do not focus further on selection in this study, noting that selective pressures are likely to be specific to a given pair or set of mtDNA types and are not generally characterised well enough to perform satisfactory inference.However, we do note that our theory gives a straightforward prediction for the functional form of mean heteroplasmy when nonzero selection is present, a sigmoid with slope set by the fitness difference (see Methods).

Probabilities of exceeding threshold heteroplasmy values
A key feature of mtDNA diseases is that pathological symptoms usually manifest when heteroplasmy in a tissue exceeds a certain threshold value, with few or no symptoms manifested below this threshold [45].The probability and timescale associated with which cellular heteroplasmy may be expected to exceed a given value is thus a quantity of key interest in clinical planning of mtDNA disease strategies.
In our model, the probability, as a function of time, of a cell containing m 1 wildtype and m 2 mutant mtDNAs can be straightforwardly derived.The resultant analytic expression involves a hypergeometric function, also an important mathematical element in expressions describing mtDNA statistics based on classical population genetics [23,46].The probability of obtaining a given heteroplasmy can therefore be computed as a sum over all copy number states that correspond to that heteroplasmy.However, as hypergeometric functions are comparatively unintuitive and computationally expensive, we here employ an approximation to the distribution of heteroplasmy based upon the above moments that are straightforwardly calculable from our model.This approximation involves fixation probabilities for each mtDNA type and a truncated Normal distribution for intermediate heteroplasmies (see Methods).In the Supplementary Information we show that this approximation corresponds well to the exact distributions calculated using the hypergeometric function.We underline that exact heteroplasmy distributions are straightfoward to compute using our approach: we use the truncated Normal approximation as it represents the exact distribution well, is more intuitively interpretable, and is computationally very inexpensive.
Using this approach, the probability with time of a cell exceeding a threshold heteroplasmy h * can be straightforwardly computed for any initial heteroplasmy, allowing rigorous quantitative elucidation of this important clinical quantity (see Methods).Fig. 4D illustrates this computation by showing the analytic probability with which thresholds h * = 0.4, 0.5, 0.6 are exceeded at a time t, given the example initial heteroplasmy h = 0.3.These results serve as a simple example of the power of our modelling approach: any other specific case can readily be addressed.Our theory thus allows general quantitative calculation of the probability (and timescale) that any given heteroplasmy threshold will be exceeded, given knowledge of the initial (or early) heteroplasmy.

Developmental sampling of embryonic heteroplasmy
We next turn to the question of estimating heteroplasmy levels in a developed organism by sampling cells during development.This principle, clinically termed preimplantation genetic diagnosis [5,47], assists in clinical planning by allowing inference of the specific heteroplasmic nature of the embryo itself rather than a population average of an affected mother's oocytes [48].However, the complicated and stochastic nature of the bottleneck makes this inference a challenging problem.
Given a heteroplasmy measurement from sampling h m , accurate preimplantation diagnosis is contingent on knowledge of the distribution P(h|h m ), that is, the probability that the embryonic heteroplasmy is h given that a measurement h m has been made.We can use our modelling framework and Bayes' theorem (see Methods) to obtain a formula for this conditional probability, allowing a rigorous probability to be assigned to inferences from preimplantation sampling.Here, as above, we employ the truncated Normal approximation for the heteroplasmy distribution, noting that the exact treatment using hypergeometric functions is straightforward but more computationally expensive.Fig. 4E illustrates this process by showing the probability distributions on embryonic heteroplasmy when measurements h m = 0.1 or 0.4 have been taken at different times during development.The increasing heteroplasmy variance through development means that substantially greater uncertainty is associated with heteroplasmy values inferred using measurements taken at later times.In conclusion, although care must be taken in applying this reasoning to cell types in which, for example, mitochondrial and cell turnover rates differ from those assumed here, or differentiation leads to tissue-specific selective factors acting on the mtDNA population, this formalism provides a general means of rigorously inferring embryonic heteroplasmy through genetic diagnosis sampling.

Discussion
We have used a general stochastic model and approximate Bayesian computation with the available experimental data on developmental mtDNA dynamics to show that the bottleneck is most likely manifest through stochastic mtDNA dynamics and partitioning, with increased random turnover later during development, a mechanism which we can describe exactly and analytically (Fig. 5).We emphasise that the bottom-up construction of our model from physical first principles both increases the flexibility and generality of our model, allowing different mechanisms to be compared together, and providing information on mtDNA dynamics throughout development rather than estimating an overall effect.We note that even though our model cannot represent the full microscopic truth underlying the mtDNA bottleneck, its ability to recapitulate the wide range of extant experimental measurements suggest that its study may yield useful insights into the effects of different treatments and perturbations on the bottleneck.
A key debate in the literature has focussed on the magnitude of the bottleneck.Some studies [11,15] have observed a depletion of mtDNA copy number during the bottleneck to minima around several hundred; other studies [18,19] have observed that mtDNA copy number remains > 10 3 .Our study shows that observed increases in heteroplasmy variance [16,12] can be achieved across this range of potential minimal mtDNA copy numbers, meaning that the much-debated magnitude of mtDNA copy number reduction is not the sole critical feature of the bottleneck, in agreement with arguments from Refs.[18,19,12].We find that the role of stochastic mtDNA dynamics can play a key role in determining heteroplasmy variance without additional mechanistic details, in keeping with approaches proposed by Ref. [11].The mechanism with the most statistical support is thus consistent with aspects from all existing proposals in the literature.
We have shown that, of the models proposed in the literature, a birth-death-partitioning model, proposed after Ref. [11] and compatible with an interpretation of Ref. [12], is the individually most likely mechanism, and capable of producing experimentally observed heteroplasmy behaviour.We cannot, given current experimental evidence, discount hybrid mechanisms, where birth-death-partitioning dominates the population dynamics but small contributions from other mechanisms provide perturbations to this behaviour, and propose experiments to conclusively distinguish between these cases (see Supplementary Information).As the expected statistics of mtDNA populations undergoing inheritance of heteroplasmic mtDNA clusters is very similar to those undergoing binomial partitioning of mtDNAs (see Supplementary Information), the inheritance of heteroplasmic nucleoids (as opposed to individual mtDNAs) is not excluded by our findings, though other recent experimental evidence suggests that this situation may be unlikely [42,41].We contend that the most likely situation may involve the partitioning of individual organelles, containing a mixture of homoplasmic nucleoids of characteristic size < 2. Notably, this case (inheritance of heteroplasmic groups, likely with fluid structure due to mixing of organellar content and mitochondrial dynamics), gives rise to statistics which our binomial model reproduces (see Supplementary Information).
As mentioned in the model description, it is likely that mitochondrial dynamics (fission and fusion of mitochondria) [35] play a role in determining natural mtDNA turnover, and particularly mtDNA turnover in the presence of pathological mutations [49], through the mechanism of mitochondrial quality control [32,37].Mitochondrial dynamics may also influence the elements of partitioning, through changes in the connectivity of the mitochondrial network.In our current model, these influences are coarse-grained into descriptions of the dynamic rates of mtDNA replication and degradation, and the characteristic elements that are partitioned at divisions.These physical parameters, as opposed to the more microscopic details of mitochondrial dynamics, are expected to be the key determinants of heteroplasmy statistics through development.Accounting for how these parameters, which summarize the relevant outputs of mitochondrial dynamics, connect to details of microscopic models of mitochondrial dynamics is an important future research direction to be addressed when more quantitative data is available.
The experimental data used to parameterise the first part of our study was taken from four studies in mice.Observation of similar dynamics in salmon [20] points towards the bottleneck being a conserved mechanism in vertebrates.We also note that our results in mice are broadly consistent with findings from recent experiments in other organisms, suggesting that in primates and humans, heteroplasmy variance may increase at early developmental stages [50,51], and that partitioning of mitochondria is binomial in HeLa cells [52].As more studies become available on human mtDNA behaviour during development we will test our model's applicability and its clinical predictions.We note that the results of a recent study of human preimplantation sampling [48] found that earlier measurements provided strong predictive power of mean heteroplasmy, about which substantial variation was recorded in the offspring -both of which results are consistent with the application of our model to theoretical sampling considerations.In addition, recent observations that the m.3243A > G mutation in humans both increases mtDNA copy number during development [53], and displays a less pronounced increase of heteroplasmy variance [51] than other mutations, are consistent with the link between heteroplasmy variance and mtDNA copy number in our theory.
The combination of modern stochastic and statistical treatments that we have employed provides a generalisable and powerful way to recapitulate experimental data and rigorously deduce underlying biological mechanisms.We have used this combination to explore pertinent questions regarding the mtDNA bottleneck (and others have used a similar philosophy to numerically explore mtDNA point mutations [25]): we hope to convince the reader that such methodology may be appropriate to explore other problems involving stochastic biological systems.We have used new experimental measurements to confirm our theoretical findings, illus-trating the beneficial and powerful coupling of mathematical and experimental approaches to address competing hypotheses in the literature.Our detailed elucidation of the bottleneck allows us to propose further experimental methodology to address the current unknowns in our theory, including the specifics of mtDNA partitioning at cell division and the roles of selective differences between mtDNA types; importantly, we also propose a strategy to investigate our claim that our most supported model is compatible with the subsetreplication picture of mtDNA dynamics.We list these experiments in full in the Supplementary Information.Finally, we believe that the theoretical foundation for mtDNA dynamics that we have produced allows increased quantitative rigour in the predictions and strategies involved in mtDNA disease therapies, illustrated by the above application of our theory to problems in mtDNA sampling strategies, disease onset timescales, and interventions to increase the power of the bottleneck.

Methods
General model for mtDNA dynamics.Our 'bottomup' model represents individual mtDNAs as elements which replicate and degrade either randomly or deterministically according to the model parameterisation.Consonant with experimental studies showing that it is often a single mutant genotype that dominates the non-wildtype mtDNA population of a cell [54], we consider two mtDNA types (wildtype and mutant), though our model can readily be extended to more mtDNA types.We denote the number of 'wild-type' mtDNAs in a cell as m 1 and the number of 'mutant' mtD-NAs as m 2 .The heteroplasmy of a cell is then h = m2 m1+m2 , that is, the population proportion of mutant mtDNA.
MtDNA dynamics within a cell cycle.Individual mtDNAs are capable of replication and degradation, with rates denoted λ and ν respectively.According to a binary categorical parameter S, these events may be deterministic (S = 0; the mtDNA population replicates and degrades by a fixed amount per unit time) or Poisson processes (S = 1; each individual mtDNA randomly replicates and degrades with average rates λ and ν).A parameter α controls the proportion of mtDNAs capable of replication: α = 1 allows all mtDNAs to replicate throughout development, α < 1 enforces a subset proportion α of replicating mtDNAs a time cutoff T after conception.
MtDNA dynamics at cell divisions.A parameter c (cluster size; a non-negative integer) dictates the partitioning of mtDNAs at cell divisions.When c = 0, partitioning is deterministic, so each daughter cell receives exactly half of its parent's mtDNA.For c > 0, partitioning is stochastic.When c = 1, partitioning is binomial: each mtDNA has a 50% chance of being inherited by either daughter cell.When c > 1, the parent cell's mtDNAs are grouped in clusters of size c before division.Each cluster is then partitioned binomially, with a 50% chance of being inherited by either daughter cell.
Different dynamic phases through development.The mtDNA population changes in different ways as develop-ment progresses, first decreasing, then recovering, then slowly growing.We include the possibility of different 'phases' of mtDNA dynamics in our model to capture this behaviour.Each phase j has its own associated pairs of λ j , ν j parameters and may either be quiescent (involving no cell divisions) or cycling (encompassing n j cell divisions).Thus, we may have an initial cycling phase with low mtDNA replication rates, so that copy number falls for several cell divisions, then a subsequent 'recovery' cycling phase with higher replication rates so that mtDNA levels are amplified, then quiescent phases as cell lineages are identified.We allow six different phases, with the first two fixed as cycling phases with the above doubling times, and the final phase fixed to include no mtDNA replication (representing the stable, final occyte state).
Initial conditions.The initial conditions of our model involve an initial mtDNA copy number m 0 (the total number of mtDNAs in the fertilised oocyte) and an initial heteroplasmy h 0 (the fraction of these mtDNAs that are mutated).
Data acquisition.We used three datasets for mtDNA copy number during mouse development: Cree [11]; Cao [18]; and Wai [12].We use two datasets for heteroplasmy variance during development: Wai [12] and Jenuth [16].By convention, we use the normalised versions of heteroplasmy variance (that is, measured variance divided by a factor h(1 − h)).Where the measurements were not given explicitly in these publications, we manually analysed the appropriate figures to extract the numerical data.For these values, we used data from correspondence regarding the Wai study (reply to Ref. [17]), and manually normalise the Jenuth dataset.The Jenuth dataset contains measurements taken in 'mature oocytes' with no time given; we assume a time of 100 dpc for these measurements, though this time is generalisable and does not qualitatively affect our results.All values are presented in the Supplementary Information.Data on cell doubling times in the mouse germ line is taken from Ref. [43], suggesting that doubling times start with an interval of every 7h, then after around 8.5 days post conception (dpc) increase to 16h, before the onset of a quiescent regime around 13.5 dpc (roughly consistent with the estimate of ∼ 25 divisions between generations in the female mouse germ line [55]).
Simulation, model selection, and parametric inference.We use Gillespie algorithms, also known as stochastic simulation algorithms [56], to explore the behaviour of our model of the bottlenecking process for a given parameterisation.For a given model parameterisation, the Gillespie algorithm is used to simulate an ensemble of 10 3 possible realisations of the time evolution of mtDNA content, and the statistics of this ensemble are recorded.The experimental data we use is derived from sets of measurements of different sizes; to compare simulation data with an experimental datapoint i corresponding to a statistic derived from n i measurements, we sampled a random subset of n i of the 10 3 simulated trajectories (all datapoints but one have n 10 3 ), and used this subset to derive the simulated statistic for comparison to datapoint i [29].
To fit the different models to experimental data we define a distance measure, a sum-of-squares residual between the E(m) (in log space) and V(h) dynamics produced by our model and observed in the data, weighted to facilitate comparison of these different quantities [29].We also constrain copy number to be < 5 × 10 5 at all points throughout development, rejecting parameterisation that disobey this criterion.Metropolis MCMC was used to identify the best-fit parameterisation according to this distance function.For statistical inference, we use approximate Bayesian computation (ABC), a statistical approach that has successfully been applied to parametric inference and model selection in dynamical systems [30] to infer posterior probability distributions both for individual models and the parameters of the models given experimental data.ABC samples posterior probability distributions on parameters that lead to behaviour within a certain threshold distance of the given data; these posteriors are shown to converge on the true posteriors as the threshold value decreases to zero (see Supplementary Information).We employed an MCMC sampler with randomly-selected initial conditions.For further details, including priors, thresholds and step sizes used in ABC, see Supplementary Information.Minimum copy number was recorded directly from the resulting trajectories; our measure of total turnover σ is defined as σ = 6 i=3 τ i ν i , the sum over quiescent dynamic phases of the product of degradation rate and phase length.
Creation of heteroplasmic mice.Heteroplasmic mice were obtained from a heteroplasmic mouse line (HB) we created previously by ooplasmic transfer [44].This mouse line contains the nuclear DNA of the C57BL/6N mouse, and mtD-NAs both of C57BL/6N and a wild-derived house mouse.Both mtDNA variants belong to the same subspecies, Mus musculus domesticus.For details on sequence divergence see Ref. [44].
Isolation and lysis of oocytes.Mice were sacrificed at the indicated ages by cervical dislocation.Ovaries were extracted and immediately placed in cryo-buffer containing 50% PBS, 25% ethylene glycol and 25% DMSO (Sigma-Aldrich, Austria) and stored at -80 o C. For oocyte extraction, ovaries were placed into a drop of cryo-buffer and disrupted using scalpel and forceps.Oocytes were collected and remaining cumulus cells were removed mechanically by repeated careful suction through glass capillaries.Prepared oocytes were then washed in PBS before they were individually placed into compartments of 96-well PCR plates (Life Technologies, Austria) containing 10 µl of oocyte-lysis buffer [50] [50 mM Tris-HCl, (p.H 8.5), 1 mM EDTA, 0.5% tween-20 (Sigma-Aldrich, Austria) and 200 µg/mL ProteinaseK (Macherey-Nagel, Germany)].Samples covered stages from primary oocytes of 3 day-old mice up to mature oocytes of 40 day-old mice.Samples were lysed at 55 o C for 2 h, and incubated at 95 o C for 10 min to inactivate Proteinase K.The cellular DNA extract was finally diluted in 190 µl Tris-EDTA buffer, pH 8.0 (Sigma-Aldrich, Austria).3µL were used per qPCR reaction.
The study was conducted according to MIQE (minimum information for publication of quantitative real-time PCR experiments) guidelines [44,60].The proportion between HB derived and C57BL/6N mtDNA was determined by ARMS-qPCR assays based on a SNP in mt-rnr2 [44].These assays were normalised to changes in the input mtDNA amount by consensus assays, located in conserved regions of mt-Co2 and mt-Co3 (see Supplementary Information).For the calculation of mtDNA heteroplasmy, the assay detecting the minor allele (C57BL/6N or wild-derived < 50%) was always used.If both specific assays gave values > 50% (which can happen around 50% heteroplasmy), the mean value of both assays was taken.All qPCR runs contained no template controls (NTCs) for all assays; these were negative in 100%.Further experimental details available in the Supplementary Information.
Analytic model.In the birth-death-partitioning model, processes within a cell cycle constitute a birth-death process which can be solved using generating functions [61].For binomial partitioning, the generating function for the system after an arbitrary number of divisions has a recursive structure [62] and an analytic solution can be obtained through solving a Riccati recurrence relation.This reasoning also extends to the different phases of replication and degradation, allowing an exact generating function to be constructed for an arbitrary point in the bottleneck.Derivatives of this generating function are then used to obtain moments of the distributions of interest.The full procedure is given in the Supplementary Information.Recall that we assume that the bottlenecking process consists of a series of dynamic phases, which may either involve cycling cells (and hence cell divisions) or quiescent cells.The expression for mean mtDNA copy number E(m, t) at time t is: where n i is the number of cell divisions in phase i (0 for quiescent phases), τ i is the length of a cell cycle in cycling phase i, τ i is the time spent in quiescent phase i (0 for cycling phases), and τ * = Σ i (n i τ i +τ i ), so that t−τ * is the time since the last cell division.E(m, t) is thus intuitively interpretable as a product of the initial copy number with the effects of halving at each cell division, and the copy number evolution through past and current cell cycles and quiescent phases.
The expression for the variance is lengthier, taking the form where Φ is a lengthy, though algebraically simple, function of all physical parameters, which we derive and present in the Supplementary Information.Once the means and variances associated with mutant and wild-type mtDNAs have been determined (for brevity, we write these as µ 1 ≡ E(m 1 , t), σ 2  1 ≡ V(m 1 , t) and µ 2 ≡ E(m 2 , t), σ 2  2 ≡ V(m 2 , t)), the relations below can be used to compute heteroplasmy statistics: Selection.The predicted mean heteroplasmy at time t assuming a constant selective pressure (though this assumption can straightforwardly be relaxed) is given by Eqn. 4, which, given Eqn.2, straightforwardly reduces to where h 0 is initial heteroplasmy and ∆λ is the increase (or decrease, if negative) in replication rate of mutant over wild-type mtDNA.Eqn.6 predicts that mean heteroplasmy in the presence of selection will follow a sigmoidal form (as expected from population dynamics [63], by the constraint that h 0 must lie between 0 and 1, and by the intuitive fact that heteroplasmy changes slow down as these limits are approached).
Threshold crossing.The probability of heteroplasmy exceeding a certain threshold h * is simply given by integrating the probability distribution of heteroplasmy between h * and 1.The exact distribution of heteroplasmy can be written as a sum over hypergeometric functions; however, for computational efficiency and interpretability, we employ an approximation to this distribution involving the truncated Normal distribution and fixation probabilities.As shown in the Supplementary Information, the distribution of heteroplasmy, taking possible fixation into account, can be well approximated by where N is the truncated Normal distribution (truncated at 0 and 1), µ and σ 2 are found numerically given our model results for E(h) and V(h), and ζ 1 ≡ P(h = 0) and ζ 2 ≡ P(h = 1) are fixation probabilities, also straightforwardly calculable from our model.The probability of threshold crossing for 0 < h * < 1 is then Inference from heteroplasmy measurements.Given a sampled measurement heteroplasmy h m , the probability P(h 0 |h m ) that embryonic heteroplasmy is h 0 is given by Bayes' theorem P(h 0 |h m ) = P(h m |h 0 )P(h 0 )/P(h m ).Assuming a uniform prior distribution on embryonic heteroplasmy (though this can be straightforwardly generalised), we thus obtain P(h 0 |h m ) = P(h m |h 0 )/ 1 0 P(h m |h 0 )dh 0 , and using the above expression for the heteroplasmy, Table 1: Source data used in this study. 1Data referenced by number of cells post-conception is assigned a time measurement assuming the 7h → 16h doubling times from Ref. [43]. 2 Mean copy number taken directly from tabulated data. 3(Weighted) average over germline cell classes presented at this time point. 4Extracted from data in figures; n not explicitly available so estimated as n = 20 from accompanying histograms and discussion. 5Manually normalised from given data. 6Data from mature oocytes in next generation: time in dpc not available. 7Extracted from data in figures in correspondence following study. 8n not explicitly available so estimated as n = 20 from accompanying histograms and discussion in original paper.
molecules (which we will refer to as 'quenched' sets) or more fluid transient colocalisations of molecules (which we refer to as 'unquenched').Furthermore, the size of these units is debated, with estimates ranging from an average size of 1.4 to 10 mtDNA molecules [39,40], and it is unknown whether the mtDNAs within a group are strictly homoplasmic or if heteroplasmic groups are possible, although current evidence, at the finest resolution, points towards homoplasmic groups of size < 2 [42,41,64,1].We will classify these different pictures with three parameters.First, the characteristic size c of an mtDNA group.Second, a classifier denoting whether these groups are quenched (in the sense that the individual constituents of a group remain the same over many cell divisions) or unquenched (in the sense that the individual constituents of a group may change between cell divisions).Third, a classifier denoting whether groups are necessarily homoplasmic, or if heteroplasmy is permitted.
An early hypothesis from Jacobs et al. [31] considered 'nucleoids' which correspond to quenched heteroplasmic groups with c > 1, retaining their internal structure across cell divisions and containing different mtDNA types.If the mitochondrial organelle is the unit of inheritance, we may expect unquenched heteroplasmic groups with c > 1, as mitochondrial dynamics act to mix the content of the mitochondrial system between cell divisions, but organelles are likely to contain more than one mtDNA molecule.If nucleoids are the units of inheritance and, as current understanding suggests, nucleoids are small and homoplasmic (if mtDNA indeed exists in groups at all), the appropriate picture is c 1, homoplasmic groups.
Here we show that the heteroplasmy statistics resulting from these different pictures of grouped inheritance collapse onto two representative cases: first, that corresponding to homoplasmic clusters with c > 1, and second, that corresponding to c = 1 (binomial inheritance).Quenching -whether mtDNA content can remix within nucleoids -is shown to be unimportant in determining heteroplasmy statistics.Our model for these different situations is as follows.We consider a cellular population as consisting of a set of mtDNA molecules, existing in groups of size c.During a cell cycle, the population of groups doubles deterministically (we ignore random birth-death dynamics in this model, in order to focus on partitioning dynamics), so that every group produces one exact copy of itself.For unquenched simulations, a new set of groups is then formed by resampling the individual mtDNA constituents of the cell.For quenched simulations (representing the situation postulated in Ref. [31]), the existing groups remain intact.At cell divisions, groups are binomially partitioned between the two daughter cells.
The model is initialised with a cell containing m 0 mtDNAs, split into (1 − h)m 0 wild-type and hm 0 mutant molecules.These mtDNAs are clustered into m 0 /c groups, according to the cluster picture under consideration (i.e.homoplasmic or heteroplasmic clusters).We simulate the subsequent doubling then partitioning of this system through cell divisions many times, assuming a constant cell cycle length, and record the cell-to-cell heteroplasmy variance with time.
Fig. 6 shows the resultant heteroplasmy variance trajectories for different cases (with h 0 = 0.1; other initial heteroplasmies showed similar behaviour).The first striking result is that the inheritance of heteroplasmic groups produces the same heteroplasmy variance as binomial partitioning, regardless of cluster size.This behaviour is due to the balance between stochasticity associated with the makeup of, and partitioning of, groups.A small number of large groups will experience substantial partitioning noise, but larger heteroplasmic groups are more likely to faithfully represent the overall cell heteroplasmy.As identified in Ref. [31], the inheritance of heteroplasmic groups thus provides a means to facilitate local mtDNA complementation while provoking no increase in heteroplasmy variance beyond that associated with binomial partitioning of elements at divisions.
We also observe that quenched populations behave in the same way as unquenched populations.In the case of homoplasmic groups, this result is obvious, as a set of homoplasmic nucleoids of a given size can only be constructed in one way for a given number of mtDNA molecules of different types.For heteroplasmic groups, this result implies that resampling the cellular population to produce a new group produces a negligible amount of additional stochasticity compared to that already present in the random makeup and inheritance of groups.Thus, the only determinant factors of heteroplasmy variance related to the inheritance of groups are whether groups are homoplasmic or heteroplasmic, and, if the former, the characteristic size of groups.
These results illustrate that the binomial inheritance model can also describe the statistics associated with heteroplasmic nucleoids of arbitrary size, over a timescale of several dozen cell divisions (suitable to describe the developmental process).The theoretical long-term behaviour of these systems involves some more subtleties.At much longer times, the probability that all mtDNA types become extinct in a cell is not negligible.When complete extinction cannot be ignored, heteroplasmy statistics become poorly defined.This extinction timescale is shorter for cluster inheritance than for binomial inheritance, as a greater variability in copy number (though not in heteroplasmy) results from each division for larger clusters.However, our simulations indicate that as long as the heteroplasmy variance associated with heteroplasmic clusters remains well defined, it matches that resulting from binomial inheritance.
We propose that a reasonable view may be that individual mitochondrial fragments, including several small, homoplasmic nucleoids, are the likely elements of inheritance at partitioning.Furthermore, there is likely some movement of these nucleoids within the mitochondrial network, and fission and fusion likely mean that a given mtDNA will not be associated with the same static mitochondrial element in perpetuity.In this case, the picture of an unquenched, heteroplasmic group of mtDNAs -those contained within a discrete element of the mitochondrial system -seems most reasonable.We can thus speculate that, as demonstrated by the previous results, the precise size of mitochondrial fragments at partitioning is not important Figure 6: Heteroplasmy variance in a model system under several different group-inheritance regimes.V(h) over many cell divisions when the elements of inheritance are heteroplasmic or homoplasmic groups of different size.Groups may be quenched (Q; constituents remain the same across cell divisions) or unquenched (UQ; constituents are randomly resampled from the cellular population each cell cycle); for homoplasmic clusters, an unquenched protocol yields identical results to the quenched protocol.V(h) behaviour differing from binomial partitioning (c = 1) is only observed for homoplasmic groups with c ≥ 2. Points for heteroplasmic groups are slightly offset in the x-direction for clarity.
for the heteroplasmy dynamics (nor indeed is whether they are quenched or unquenched).Our simple binomial partitioning model is thus consistent with what one might consider the most physiologically plausible model, and indeed with any models not involving large and strictly homoplasmic groups as the elements of mitochondrial inheritance.

Parametric inference for bottlenecking dynamics
Our model is a function of the parameter set1 θ = {ν i , λ i , n i , τ i , S, α, T, c, h 0 , m 0 , δλ}.For the following parameters we use uninformative uniform priors on the given interval: 10 6 ].The following values are fixed from experimental studies [43,55]: n 1 = 29; n 2 = 7; τ 1 = 7hr; τ 2 = 16hr.λ 6 = 0hr −1 is fixed to avoid mtDNA proliferation after development; h 0 = 0.2 is fixed as an intermediate value as heteroplasmy variance measurements are generally normalised; δλ, a parameter allowing a difference in replicative rates between mutant and wildtype mtDNAs, is fixed at zero throughout as we ignore selective pressure.The parameter τ i for i > 2 is used to determine the length of time spent in different quiescent phases and is subject to the uniform prior τ i ∈ [0, 50]day.
Given these priors, we use an approximate Bayesian computation (ABC) approach to build a posterior distribution over the parameters in our bottlenecking model [30].ABC involves using a summary statistic ρ(θ, D) to compare the available data D to the predictions of a model given parameters θ.If parameter sets are sampled from the set for which ρ ≤ , where is a threshold difference between the resulting model behaviour and experimental data, the posterior distribution P (θ|ρ(θ, D ≤ )) is sampled, which is argued to sufficiently approximate P (θ|D) for suitably small [65].
We define a residual sum-of-squares difference between the results of a simulated model and experimental data [29]: where D denotes experimental data.We thus amalgamate experimental results of two types: mean mtDNA copy number (with N m data points measuring E D (m) at times t (i) D,m ), and mean and variance of heteroplasmy (with N h data points measuring V D (h) at times t (i) D,h ).The sets of data for E(m) and V(h) contain different numbers of points and are of different absolute magnitudes.We compensate for these differences by using the logarithms of copy number measurements (as these values span several orders of magnitude), and weighting parameter A 1 = 10 3 .This weighting parameter compensates for the different magnitudes and number of datapoints in each class of measurement, ensuring that the contribution to the total residual from each set of data is of comparable magnitude.Our summary statistic thus records a residual sum-of-squares difference between experiment and simulation values for log E(m) and V(h) at each time point where an experimental measurement exists.We performed our model selection process using several different alternative protocols, including comparing logarithms of V(h) measurements (in contrast to the raw values) and varying A 1 over orders of magnitude from 10 2 − 10 4 (corresponding to unbalanced weighting, favouring E(m) and V(h) data respectively).In all cases, the BDP model identified in the Main Text experienced substantially more support than any alternative.For inference involving the new dataset from the HB model system, we use the default protocol above and set A 1 = 3 × 10 3 to account for the threefold decrease in available V (h) datapoints.
We use an MCMC implementation of ABC, whereby we construct a Markov chain θ i , where each state consists of a set of trial parameters to be assessed.We create θ i+1 by perturbing each parameter within θ i with a perturbation kernel consisting of a Normal distribution on each parameter with standard deviations between 0.1 − 1% of the width of the prior (varied as the model depends more sensitively on some parameters than others).In the case of discrete parameter c, a continuous representation c is used and varied in the MCMC approach, with c = 100c .We accept θ i+1 as the new state of the chain if ρ(θ i+1 , D) ≤ .We ran 10 6 MCMC iterations for ABC model selections and checked convergence by running five instances of each simulation for different random number seeds.
For the initial optimisation of model fitting, we ran 10 6 MCMC steps using the protocol above but accepting a move according to the Metropolis-Hastings protocol [66], recording the parameterisation leading to the lowest recorded residual.In this case we used uninformative initial conditions, with identical choices for all rate parameters, corresponding to an inaccurate trajectory of copy number and heteroplasmy variance.For model selection, we used the protocol above, with a different set of parameters θ M for each model M , with each MCMC step proposing a random model from the Cao alone, Wai alone, and BDP set described in the text, as in the SMC ABC model selection protocol proposed in Ref. [30].We record the proportion of accepted steps involving each model type.The parameterisations found through initial optimisation were used as initial conditions in the ABC model selection and inference simulations.
Initial optimisation identified parameterisations all displaying residuals under = 50.We chose = {45, 50, 60, 75} for the ABC model selection simulations to display the varying degrees of support for each model as stricter agreement with experiment was enforced.We chose = 40 for the ABC inference of BDP model parameterisation to ensure these models all displayed better fits to data than the alternative models.In Fig. 7 we illustrate the distribution of squared residuals for the BDP model under a range of values.

Posteriors for all variables and datasets
In Fig. 8 we display all posterior distributions for all parameters resulting from our ABC approach assuming the BDP model.There is substantial variability in the possible timescales and magnitudes of random turnover associated with each random dynamic phase i > 2, exemplified by the complicated and bimodal structure of the posteriors on these parameters.This variability reflects the fact that an increase in heteroplasmy variance can be achieved through a variety of specific mtDNA trajectories, and current experimental data is insufficient to distinguish specific time behaviours within this variety.However, the total contribution of each random phase to the overall dynamics is more constrained, as shown in the posterior distribution on a measure of total random turnover σ = 6 i=3 τ i ν i .This quantity is the sum over all later phases of the product of the length of that phase and the rate of random turnover, thus giving a measure of total random turnover.The fact that this posterior is more tightly constrained than the posteriors on individual t i , ν i parameters suggests that the required mtDNA turnover can be achieved through a range of specific dynamic trajectories from the inferred mechanism: for example, the exact time at which random mtDNA turnover sharply increases is currently flexible (though constrained to lie around 25 dpc) without more detailed data.This flexibility is also observed in the trajectories of posterior distributions in the main text.

Experimental measurements
Table 2 contains the measurements of heteroplasmy h, mean heteroplasmy E(h), raw and normalised heteroplasmy variance V(h) and V (h), and number of datapoints n, from the HB model system.Experimental procedures are described in Methods; further specifics follow.
Master-mixes for triplicate qPCR reactions contained 1x buffer B (Solis BioDyne, Estonia); 4.5 MgCl2 for the ARMS and the Co3 consensus assays, and 3.5 mM MgCl2 for the Co2 consensus assay; 200 M of each of the four deoxynucleotides (dNTPs, Solis BioDyne, Estonia), HOT FIREPol DNA polymerase according to the manufacturers instructions (Solis BioDyne, Estonia), 300 nM of each primer and 100 nM hydroloysis probe (Sigma-Aldrich, Austria).Per reaction 12 µL of master-mix and 3 µL DNA were transferred in triplicates to 384-well PCR plates (Life Technologies, Austria) using the automated pipetting system epMotion 5075TMX (Eppendorf, Germany).Amplification was performed on the ViiA 7 Real-Time PCR System using the ViiA T M 7 Software v1.The standard curve method was applied.Amplification efficiencies were determined for each run separately by DNA dilution series consisting of DNA from wild-derived mice, harbouring the respective analysed mtDNA.Typical results: slope = -3.665,-3.562, -3.461, -3.576; mean efficiency = 0.87, 0.9, 0.94, 0.90; and Y-intercept = 32.2,28.3, 33.8, 34.5; for the consensus Co2, consensus Co3, C57BL/6N and HB assays respectively.Coefficient of correlation was > 0.99 in all assays in all runs.All target samples lay within the linear interval of the standard curves.To test for specificity, in each run a negative control sample, i.e. a DNA sample of a mouse harbouring the mtDNA of the non-analysed type in the heteroplasmic mouse (i.e.C57BL/6N or HB mtDNA) was measured.All assays could discriminate between C57BL/6N and HB mtDNA at a minimum level of 0.2%.Target sample DNA was tested for inhibition by dilution in Tris-EDTA buffer (Sigma-Aldrich, Austria), pH 8.0.

Bottlenecking mechanisms and further experimental elucidation
Here we summarise potential mechanisms for the bottleneck that conflict with our statistical interpretation, highlighting the reasons for the conflict.We also propose further experiments that would efficiently provide more evidence to distinguish these hypotheses.

Other proposed mechanisms
Random partitioning of homoplasmic mtDNA clusters.Ref. [18] suggests a less powerful depletion of mtDNA copy number during early development than assumed by other studies, with heteroplasmy variance increase instead being explained by the partitioning of clusters of mtDNA at cell divisions.However, the time period over which Refs. [12] and [16] observe increasing heteroplasmy variance corresponds to a situation in which germ line cells are largely quiescent, immediately suggesting that partitioning at cell divisions cannot explain increasing variance (as cell divisions do not occur).Furthermore, results from our model suggest that, unless these clusters are very small, this mechanism would immediately lead to a rather higher and sharper increase in heteroplasmy variance than observed.
Replication of a specific subset of mtDNAs during folliculogenesis.Ref. [12] proposes a mechanism in which only a subset of mtDNAs replicate during folliculogenesis.There are several specific dynamic schemes by which this mechanism could be manifest.The first that we consider involves the following scenario: at some point during development, around the start of folliculogenesis, a specific subset of mtDNAs in each cell is 'marked' as able to replicate (we next consider the case in which this subset is more plastic with time).In this case, the effect of 'switching off' replication of a subset of mtDNAs depends on the balance of replication and degradation rates of the mtDNA population: • Low replication, low degradation.In this case, the population stays largely static; the switching off of replication has little effect, and the heteroplasmy variance cannot increase to the levels observed in experiment.
• Low replication, high degradation.In this case, the high degradation rate ensures that the non-replicating mtDNAs are removed from the cell, providing a 'bottleneck' as only the replicating mtDNAs remain.However, this regime yields a transient period of mtDNA copy number depletion, while the non-replicating mtDNAs are degrading but the (small) population of replicating agents remains low.This copy number depletion is not observed.
• High replication, high degradation.In this case, non-replicating mtDNAs are removed and replicating mtDNAs are capable of fast enough replication to survive the transient drop in copy number.However, the rates associated with this mechanism are necessarily high enough such that the increase in heteroplasmy variance is very sharp, notably more so than the smooth increase with time observed in experiment (see, for example, Fig. 2A in the Main Text).
Measurement Purpose MtDNA copy number before and after cell divisions and/or variance of copy number between daughter cells To elucidate mechanism of mtDNA partitioning and whether this partitioning is deterministic or stochastic.

Copy number trajectories with different mtDNA heteroplasmies
To assess the modulation of copy number dynamics by mtDNA heteroplasmy via retrograde signalling.Measurement of mean heteroplasmy through development, with a variety of mtDNA type pairings

To assess and quantify to what extent selection modulates mtDNA dynamics during germline development Copy number measurements after upregulation of mitophagy
To assess the presence and strength of compensatory mechanisms that may act to preserve mtDNA copy numberand hence whether upregulating mitophagy will act to increase mtDNA turnover or simply lower copy number.

Heteroplasmy variance after upregulation of mitophagy
To assess the efficacy of mitophagy for increasing the power of the bottleneck.Heteroplasmy distribution in cells after the bottleneck from sampled/known initial heteroplasmy To confirm predictions for threshold crossing and statistics between generations.

BrU incorporation in oocytes between 30 and 40 dpc
To confirm the random turnover mechanism: we expect a large proportion of BrU incorporation subset of mtDNAs to be observed in this time period (see Section 6.2).

Mitochondrial ultrastructure and mtDNA localisation during development
To assess and characterise any potential modulation of the size of units of mitochondrial inheritance by mitochondrial dynamics through development, in particular, investigating whether there is time-varying modulation of cluster size at points of division.
Table 3: Experiments for further elucidation of the mtDNA bottleneck.

Experimental elucidation
In Table 3 we list several classes of potential experimental protocols that would assist in further elucidation of the bottlenecking mechanism and our predictions.Potentially useful results include further characterisation of the microscopic detail underlying mtDNA dynamics during development, confirmation of our random turnover model, assessing degree to which heteroplasmy modulates copy number dynamics and exploring our predictions relating mitophagy and bottlenecking power.

Mitophagy regulation
The results from our model suggest a potential clinical pathway for increasing heteroplasmy variance, and thus the power of the bottleneck to remove heteroplasmic cells.We have shown that upregulation of mtDNA degradation (for example, through increasing mitophagy) leads to lower mtDNA copy numbers and greater heteroplasmy variance.It is unclear whether a given treatment will have the sole effect of upregulating mitophagy: it seems likely that compensatory mechanisms (which we do not explicitly model, but may include retrograde signalling [67]) will engage to stabilise mtDNA copy number.However, such mechanisms would most straightforwardly be expected to act through increasing mtDNA proliferation, thus having the net effect of increasing mtDNA turnover.We have shown that such an increase in turnover also increases the heteroplasmy variance in a population.We therefore propose that upregulating mitophagy may be a fruitful pathway of investigation for increasing bottlenecking power, either as a standalone effect or due to the action of compensatory mechanisms it may invoke.
Speculatively, potential strategies to upregulate mitophagy may include the limited use of uncouplers to accelerate the mitophagy normally involved in quality control [68]; targetted chemical treatments with agents that have been identified as regulating mitophagy, including glutathione in yeast [69] and C 18 -pyridium ceramide in human cancer cells [70]; modulation of mitochondrial ultrastructure and dynamics to upregulate fission, intrinsically linked to the process of mitophagy [37,71]; or the use of existing drugs which have been found to modulate mitophagy, such as Efavirenz [72].

Heteroplasmy statistics
We have defined heteroplasmy by To find statistics for this quantity we consider the Taylor expansion of a function f (X 1 , X 2 ) of two random variables X 1 , X 2 about a point (µ 1 , µ 2 ), where µ i = E(X i ).We assume that the moments of X i are well-defined and both have zero probability mass at X i = 0.The Taylor expansion is: G 0 (z, t|m 0 ) = (z − 1)νe (λ−ν)t − λz + ν (z − 1)λe (λ−ν)t − λz + ν m0 (30) where the 0 subscript signifies that no divisions have occurred, and we have specifically labelled the base of G 0 as g 0 for later convenience.
Generating function over cell divisions.We now consider a system undergoing cell divisions.Now, we have a population of organelles with time evolution described by a generating function G = [g] m0 and subject to binomial partitioning at cell division.The probability distribution of m after a single cell division is: where m i,a , m i,b mean respectively the number of individuals after and before the ith cell division, and the subscript in P 0 denotes the fact that this function refers to time evolution within a cell cycle (with no division).The sum takes into account all possible configurations of the system up to the cell division then all possible configurations afterwards, with weighting according to a binomial partitioning.This line of reasoning can straightforwardly be extended to n cell divisions [62]: where Φ i is a 'probability propagator' of the form and m n+1,a ≡ m 0 .For clarity, we introduce the nomenclature: where we have used the identity .Comparing Eqns.38 and 41 and following this process by induction we can see that the overall generating function is G n = h m0 0 , where h is the solution to the recursive system We note that Φ is just a notational simplification and is straightforwardly calculable by inserting Eqns.56-59 into Eqn.64 then computing Eqn.73.
Constant population size.For generality, we consider enforcing a constant population size in post-mitotic cells (not undergoing divisions).This process involves setting λ = ν, so the net gain in mtDNA is zero.If we write λ = ν + and take the limit → 0, Eqn. 30 becomes To enforce a constant mean population size in mitotic cells, it is necessary to balance the expected loss of mtDNA through repeated divisions with an expected increase during the cell cycle.This balance can be accomplished by setting λ = ν + ln 2 τ .Writing λ = ν + ln 2 τ + and taking the → 0 limit we obtain In both these cases, the same approach as above can be used to derive moments of the resulting probability distributions.Explicit distributions.The probability of observing exactly m mtDNAs of a given type can be found from the generating function with The distribution of heteroplasmy is then given by where I(h , h) is an indicator function returning 1 if h = h and 0 otherwise.Computing the probability of observing a given heteroplasmy thus involves a sum, over all mtDNA states that correspond to that heteroplasmy, of the probability of that state.
The evaluation of hypergeometric functions is more computationally demanding than that of more common mathematical functions, and the infinite sums at first glance seem intractable.However, in practise and using parameterisations from our inferential approach, vanishingly little probability density exists at m 1 , m 2 > 5 × 10 5 , corresponding to the biological observation that mtDNA copy number is very unlikely to exceed this value.Dynamic programming then allows these sums to be performed straightforwardly.
Finally, the computation of P(m = 0, t) is important in our analysis of the characterisation of key distributions using the first two moments (see below), where it appears as P(m 2 = 0, t), the probability of wildtype fixation.This is relatively straightforward to address analytically as when m = 0, Eqn.77 reduces to P(0, t) = G overall | z=0 , which in the notation above is simply:

Figure 1 :
Figure 1: The mitochondrial bottleneck, and elements of a general model for bottlenecking mechanisms.(A) The mtDNA bottleneck acts to produce a population of oocytes with varying heteroplasmies from a single initial oocyte with a specific heteroplasmy value.During development, mtDNA copy number per cell decreases (by a debated amount, which we address; see Main Text) then recovers, suggesting a 'bottleneck' of cellular mtDNA populations.(B) Cellular mtDNA populations during the bottleneck are modelled as containing wildtype and mutant mtDNAs.MtDNAs can replicate and degrade within a cell cycle, with rates λ and ν respectively.(C) At cell divisions, the mtDNA population is partitioned between two daughter cells either deterministically, binomially, or through the binomial partitioning of mtDNA clusters.(D) Symbols used to represent quantities and model parameters used in the main text, and their biological interpretations.

Figure 2 :
Figure 2: Different mechanisms for the mtDNA bottleneck.

Figure 3 :
Figure 3: Parameterisation of the BDP model and inferred details of bottleneck mechanism.Trajectories of (A) mean copy number E(m) and (B) normalised heteroplasmy variance V (h) resulting from BDP model parameterisations sampled using ABC with a threshold = 40.* denotes measurements in mature oocytes, modelled as 100 dpc (see Methods).Note: the range in (B) does not correspond to a credibility interval on individual measurements, but rather on an expected underlying (population) variance, from which individual variance measurements are sampled.We thus expect to see, for example, several measurements lower than this range due to sampling limitations (see text).(C) Posterior distributions on mtDNA turnover ν with time.(D) Posterior distribution on min E(m), the minimum mtDNA copy number reached during development.(E) Posterior distribution on σ = 6i=3 τ i ν i , a measure of the total amount of mtDNA turnover.
Fig. A and the mean curve in Fig. C arises because Fig.A shows the behaviour of the model under a single, optimised parameterisation, whereas Fig. C

C
HB theory: 95% c.i. HB theory: s.d.HB theory: mean HB experiments Prev (BALB/c) theory: 95% c.i. Predictions and experimental verification of the BDP model.(A) New V (h) measurements from the HB mouse system, with optimised fits for the BDP, Wai (b) and Cao models.(B) Posterior probabilities of each model given this data under decreasing ABC threshold: = {50, 40, 30, 25}.(C) All V (h) measurements from the HB model (points) with inferred V (h) behaviour from ABC applied to the BDP model (red curves).As in Fig. , this range does not correspond to a credibility interval on individual measurements, but rather on an expected underlying (population) variance, from which individual variance measurements are sampled.The inferred behaviour strongly overlaps with the inferred behaviour for the BALB/c system (blue curves), suggesting that the BDP model applies to a genetically diverse range of systems.(D) Heteroplasmy distributions.The transformation

Figure 4 :
Figure 4: Quantitative influences and clinical results from our bottlenecking model.(A-C) Trajectories of copy number E(m) and normalised heteroplasmy variance V (h) resulting from perturbing different physical parameters.Trajectory C labels the 'control' trajectory resulting from a fixed parameterisation; black dots show experimental data; * denotes measurements from primary oocytes, modelled at 100 dpc.(A) Increasing (T + ) and decreasing (T − ) mtDNA turnover (both mtDNA replication and degradation) by 20%.(B) Increasing (M + ) and decreasing (M − ) mtDNA degradation throughout development by a constant value (2 × 10 −4 , in units of day −1 ), while keeping replication constant.(C) Applying a positive (S + ) and negative (S − ) selective pressure to mutant mtDNA by 5 × 10 −6 day −1 .(D) Probability of crossing different heteroplasmy thresholds h * with time, starting with initial heteroplasmy h 0 = 0.3.(E) Probability distributions over embryonic heteroplasmy h 0 given a measurement hm from preimplantation sampling (** hm = 0.1; *** hm = 0.4) at different times.

Figure 5 :
Figure5: Model for the mtDNA bottleneck.A summary of our findings.(A) There is most statistical support for a bottlenecking mechanism whereby mtDNA dynamics is stochastic within a cell cycle, involving random replication and degradation of mtDNA, and mtDNAs are binomially partitioned at cell divisions.(B) This mechanism results in heteroplasmy variance increasing both due to stochastic partitioning at divisions and due to random turnover.The absolute magnitude of the copy number bottleneck is not critical: a range of bottleneck sizes can give rise to observed dynamics.Random turnover of mtDNA increases heteroplasmy variance through folliculogenesis and germline development.

Figure 7 :
Figure 7: Residual distributions at different ABC thresholds .The distribution of squared residuals corresponding to individual experimental datapoints compared to an ensemble of simulated trajectories for (top) log E(m) (bottom) V(h).The V(h) residuals are scaled by A 1 = 10 3 to ensure that the two sets of measurements are compared on a quantitatively equal footing.As is decreased ( 1,2,3,4 = 40, 50, 75, 100), distributions of residuals from accepted trajectories tighten around zero.

1 (
Life Technologies, USA).DNA denaturation and enzyme activation were performed for 15 min at 95 o C. DNA was amplified over 40 cycles consisting of 95 o C for 20 sec, 58 o C for 20 sec and 72 o C for 40 sec for both assays.

Figure 8 :
Figure 8: Posterior distributions on model parameters.The posterior distributions on individual model parameters, assuming the inferred BDP bottlenecking mechanism.Replication rates are presented as κ = λ − ν, thus representing overall proliferation rates of mtDNA.Units are omitted for clarity.Pale, single-values distributions correspond to parameter values fixed within the model (κ 6 = 0 to prevent mtDNA proliferation after development; τ 1 = 7hr, τ 2 = 16hr fixed by data on cell doubling times; h 0 = 0.2 fixed for simplicity as heteroplasmy variances are normalised; δλ = 0 fixed to avoid varying selective pressure).The 'turnover' parameter, described in the text, is

Table 2 :
New heteroplasmy measurements from the HB model system.Heteroplasmy measurements and statistics from the HB model system.Ages are given in days after birth.