Mitochondrial network fragmentation modulates mutant mtDNA accumulation independently of absolute fission-fusion rates

Mitochondrial DNA (mtDNA) mutations cause severe congenital diseases but may also be associated with healthy aging. MtDNA is stochastically replicated and degraded, and exists within organelles which undergo dynamic fusion and fission. The role of the resulting mitochondrial networks in determining the time evolution of the cellular proportion of mutated mtDNA molecules (heteroplasmy), and cell-to-cell variability in heteroplasmy (heteroplasmy variance), remains incompletely understood. Heteroplasmy variance is particularly important since it modulates the number of pathological cells in a tissue. Here, we provide the first widereaching mathematical treatment which bridges mitochondrial network and genetic states. We show that, for a range of models, the rate of increase in heteroplasmy variance, and the rate of de novo mutation, is proportionately modulated by the fraction of unfused mitochondria, independently of the absolute fissionfusion rate. In the context of selective fusion, we show that intermediate fusion/fission ratios are optimal for the clearance of mtDNA mutants. Our findings imply that modulating network state, mitophagy rate and copy number to slow down heteroplasmy dynamics when mean heteroplasmy is low, could have therapeutic advantages for mitochondrial disease and healthy aging.


Introduction
Mitochondrial DNA (mtDNA) encodes elements of the respiratory system vital for cellular function.Mutation of mtDNA is a leading hypothesis for the cause of normal aging (López-Otín et al., 2013;Kauppila et al., 2017), as well as underlying a number of heritable mtDNA-related diseases (Schon et al., 2012).Cells typically contain hundreds, or thousands, of copies of mtDNA per cell: each molecule encodes crucial components of the electron transport chain, which generates energy for the cell in the form of ATP.Consequently, the mitochondrial phenotype of a single cell is determined, in part, by its fluctuating population of mtDNA molecules.The broad biomedical implications of mitochondrial DNA mutation, combined with the countable nature of mtDNAs and the stochastic nature of their dynamics, offer the opportunity for mathematical understanding to provide important insights into human health and disease.
An important observation in mitochondrial physiology is the threshold effect, whereby cells may often tolerate relatively high levels of mtDNA mutation, until the fraction of mutated mtDNAs (termed heteroplasmy) exceeds a certain critical value where a pathological phenotype occurs (Rossignol et al., 2003;Picard et al., 2014;Stewart and Chinnery, 2015;Aryaman et al., 2017).Fluctuations within individual cells mean that the fraction of mutant mtDNAs per cell is not constant within a tissue (Figure 1A), but follows a probability distribution which changes with time (Figure 1B).In the following, motivated by a general picture of aging, we will largely focus on the setting of non-dividing cells, which possess two mtDNA variants (we will consider de novo mutation in Box 1).The variance of the distribution of heteroplasmies informs upon the fraction of cells Figure 1.A simple model bridging mitochondrial networks and genetics yields a wide-reaching, analytically obtained, description of heteroplasmy variance dynamics.(A) A population of cells from a tissue exhibit inter-cellular heterogeneity in mitochondrial content: both mutant load (heteroplasmy) and copy number.(B) Inter-cellular heterogeneity implies that heteroplasmy is described by a probability distribution.Cells above a threshold heteroplasmy (h * , black dashed line) are thought to exhibit a pathological phenotype.The low-variance distribution (black line) has fewer cells above a pathological threshold heteroplasmy than the high-variance distribution (red line).(C) The chemical reaction network we use to model the dynamics of mitochondrial DNA (see Main Text for a detailed description).MtDNAs are assigned a genetic state: mutant (M ) or wild-type (W ), and a network state: singleton (i.e.unfused, S) or fused (F ).(D) The central result of our work is, assuming that a cell at time t = 0 is at its (deterministic) steady-state, heteroplasmy variance (V(h)) approximately increases with time (t), mitophagy rate (µ) and the fraction of mitochondria that are unfused (fs), and decreases with mtDNA copy number (n).Importantly, V(h) does not depend on the absolute magnitude of the fission-fusion rates.
above a given pathological threshold (Figure 1B).Therefore heteroplasmy variance is related to the number of dysfunctional cells within a tissue, and is a quantity of importance in addition to mean heteroplasmy.
Mitochondria exist within a network which dynamically fuses and fragments.Although the function of mitochondrial networks remains an open question (Hoitzing et al., 2015), it is often thought that a combination of network dynamics and mitochondrial autophagy (termed mitophagy) perform quality control on the mitochondrial population (Twig et al., 2008).However, the observation of pervasive intra-mitochondrial mtDNA mutation (Morris et al., 2017) and universal heteroplasmy in the human population (Payne et al., 2012) suggest that mitochondrial quality control mechanisms at the genetic level are weak, at least for some mitochondrial mutations.Neutral genetic models are therefore a rational means to understand mtDNA dynamics, whereby each species has the same probability per unit time of degradation, and neither species has a replicative advantage (Chinnery and Samuels, 1999;Poovathingal et al., 2009;Johnston and Jones, 2016).
A number of studies have attempted to understand the impact of the mitochondrial network on mitochondrial dysfunction through computer simulation (reviewed in (Kowald and Klipp, 2014)).Mouli et al. (2009) use a stochastic model to investigate the role of the mitochondrial network in quality control, where mitochondria incur "damage" at a constant rate, which was linked sigmoidally to a qualitative mitochondrial "activity".The authors investigated selective fusion (where fusion occurs between sufficiently active mitochondria), and assumed selective mitophagy (where only mitochondria with sufficiently low activity are degraded).The authors concluded that i) a higher fusion rate means more effective clearance of damaged mitochondria when fusion is selective; ii) intermediate fusion rates are optimal for selective mitophagy.Patel et al. (2013) take a similar approach, simulating autophagy, replication, fusion, fission and transport of mitochondria, finding that when network dynamics are coupled to transport, that transport can indirectly modulate mitochondrial health through mitochondrial dynamics.In another pair of studies by Tam et al. (2013Tam et al. ( , 2015) ) the authors applied a chemical reaction model where mutant and wild-type mtDNAs undergoing selective fusion, fission, mitophagy and replication, where the cell is partitioned into spatial compartments.The authors argue that while more frequent fission-fusion events increase the occurrence of mutant-rich mitochondria, which may be removed by selective mitophagy, these mitochondria exist for a shorter period of time, hence the existence of a trade-off (Tam et al., 2015).The authors also argue that slower fission-fusion rates result in slow spatial mixing of mtDNA, implying greater stochasticity, and a greater propensity for the expansion of mutant mtDNAs (Tam et al., 2013).Finally, a study by Figge et al. (2012) which models transitions within an abstract "quality" space, argues that if mitochondrial fission induces damage, decelerating fission-fusion cycles may improve mitochondrial quality; note that mitochondrial fission is usually interpreted as exposing damaged mitochondria rather than itself being harmful (Twig et al., 2008;Kowald and Klipp, 2014).
Previous attempts to link mitochondrial genetics and network dynamics, while important for breaking ground, have centred around complex computer simulations, making it difficult to deduce general laws and principles.In contrast, we take a simpler approach in terms of our model structure (Figure 1C), allowing us to derive explicit, interpretable, mathematical formulae which provide intuitive understanding, and give a direct account for the phenomena which are observed in our model (Figure 1D).Our results hold for a range of variant model structures.Simplified approaches using stochastic modelling have shown success in understanding mitochondrial physiology from a purely genetic perspective (Chinnery and Samuels, 1999;Capps et al., 2003;Johnston and Jones, 2016).Furthermore, there currently exists limited evidence for pronounced, universal, selective differences of mitochondrial variants in vivo (Stewart and Larsson, 2014;Hoitzing, 2017).Our basic approach therefore also differs from previous modelling attempts, since our model is neutral with respect to genetics (no replicative advantage or selective mitophagy) and the mitochondrial network (no selective fusion).Evidence for negative selection of particular mtDNA mutations has been observed in vivo (Ye et al., 2014;Morris et al., 2017); we therefore extend our analysis to explore selectivity in the context of mitochondrial quality control using our simplified framework.
Here, we reveal the first general mathematical principle linking network state and heteroplasmy statistics (Figure 1D).This link shows analytically, for a broad range of situations, that the expansion of mtDNA mutants is strongly modulated by network state, such that the rate of increase of heteroplasmy variance, and the rate of accumulation of de novo mutation, is proportional to the fraction of unfused mitochondria.This result stems from the notion that fusion shields mtDNAs from turnover, since autophagy of large fragments of the mitochondrial network are unlikely, which effectively rescales time.Importantly, we show that heteroplasmy variance is independent of the absolute magnitude of the fusion and fission rates due to a separation of timescales between genetic and network processes (in contrast to (Tam et al., 2015)).Furthermore, we find that when fusion is selective, intermediate fusion/fission ratios are optimal for the clearance of mutated mtDNAs (in contrast to (Mouli et al., 2009)).When mitophagy is selective, complete fragmentation of the network results in the most effective elimination of mitochondrial mutants (in contrast to (Mouli et al., 2009)).We also confirm that mitophagy and mitochondrial DNA copy number also affect the rate of accumulation of de novo mutations (Johnston and Jones, 2016).We suggest that pharmacological interventions which promote fusion, slow mitophagy and increase copy number during youth may slow the rate of accumulation of pathologically mutated cells, with implications for mitochondrial disease and aging.

Results
Stochastic modelling of the coupling between genetic and network dynamics of mtDNA populations Our modelling approach takes a chemical master equation perspective by combining a general model of neutral genetic drift (for instance, see (Chinnery and Samuels, 1999;Johnston and Jones, 2016)) with a model of mitochondrial network dynamics.We consider the existence of two mitochondrial alleles, wild-type (W ) and mutant (M ).MtDNAs exist within mitochondria, which undergo fusion and fission.We therefore assign mtDNAs a network state: fused (F ) or unfused (we term "singleton", S).This representation of the mitochondrial network allows us to include the effects of the mitochondrial network in a simple way, without the need to resort to a spatial model or consider the precise network structure, allowing us to make analytic progress and derive interpretable formulae in a more general range of situations.
Our model can be decomposed into three notional blocks (Figure 1C).Firstly, the principal network processes denote fusion and fission of mitochondria containing mtDNAs of the same allele 1) 3) where X denotes either a wild-type (W ) or a mutant (M ) mtDNA (therefore a set of chemical reactions analogous to (Equation 1)-(Equation 3) exist for both DNA species).γ and β are the stochastic rate constants for fusion and fission respectively.Note that we constrain all fusion reactions to have the same rate constant, since this is the simplest model of fusion.Secondly, mtDNAs are replicated and degraded through a set of reactions termed genetic processes.A central assumption is that all degradation of mtDNAs occur through mitophagy, and that only small pieces of the mitochondrial network are susceptible to mitophagy; for parsimony we take the limit of only the singletons being susceptible to mitophagy 5) 6) where λ and µ are the replication and mitophagy rates respectively, which are shared by both W and M resulting in a so-called 'neutral' genetic model.∅ denotes removal of the species from the system.The effect of allowing non-zero degradation of fused species is discussed in STAR Methods (see (Equation 82) and Figure S2G).Note that replication of a singleton changes the network state of the mtDNA into a fused species, since replication occurs within the same membrane-bound organelle.An alternative model of singletons which replicate into singletons leaves our central result (Figure 1D) unchanged (see (Equation 81)).The system may be considered neutral since both W and M possess the same replication and degradation rates per molecule of mtDNA at any instance in time.
Finally, mtDNAs of different genotypes may interact through fusion via a set of reactions we term network cross-processes: 7) 8) 9) We note that any fusion or fission event which does not involve the generation or removal of a singleton leaves our system unchanged; we term such events as non-identity-changing processes, which can be ignored in our system (see STAR Methods for a discussion of rate renormalization).Note that we have neglected de novo mutation in the model description above, see Box 1 for a treatment of de novo mutation using a modified infinite sites Moran model.We found that treating λ = const led to instability in total copy number (see STAR Methods), which is not credible.We therefore favoured a state-dependent replication rate such that copy number is controlled to a particular value, as has been done by previous authors (Chinnery and Samuels, 1999;Capps et al., 2003;Johnston and Jones, 2016).Allowing lower-case variables to denote the copy number of their respective molecular species, we will focus on a linear replication rate of the form (Hoitzing et al., 2017;Hoitzing, 2017): 10) where w T = w s + w f is the total wild-type copy number, and similarly for m T .b is a parameter which determines the strength with which total copy number is controlled to a target copy number, and κ is a parameter which is indicative of (but not equivalent to) the steady state copy number.δ indicates the relative contribution of mutant mtDNAs to the control strength and is linked to the "maintenance of wild-type" hypothesis (Durham et al., 2007;Stewart and Chinnery, 2015).When 0 ≤ δ < 1, and both mutant and wild-type species are present, mutants have a lower contribution to the birth rate than wild-types.When wild-types are absent, the population size will be larger than when there are no mutants: hence mutants have a higher carrying capacity in this regime.We have modelled the mitophagy rate as constant per mtDNA.We do, however, explore relaxing this constraint below by allowing mitophagy to be a function of state, and also affect mutants differentially under quality control.Analogues of this model (without a network) have been applied to mitochondrial systems (Chinnery and Samuels, 1999;Capps et al., 2003).Overall, our simple model consists of 4 species (W S , W F , M S , M F ), 6 parameters and 15 reactions, and captures the central property that mitochondria fragment before degradation (Twig et al., 2008).

Mitochondrial network state rescales the linear increase of heteroplasmy variance over time independently of fission-fusion rate magnitudes
We first performed a deterministic analysis of the system presented in (Equation 1)-(Equation 10), by converting the reactions into an analogous set of four coupled ordinary differential equations (see (Equation 43)-(Equation 46)), and choosing a biologically-motivated approximate parametrization (which we will term as the 'nominal' parametrization, see STAR Methods).Figures 2A-B show that copy numbers of each individual species change in time such that the state approaches a line of steady states ((Equation 48)-(Equation 50)), as seen in other neutral genetic models (Capps et al., 2003;Hoitzing, 2017).Upon reaching this line, total copy number remains constant (Figure S2A) and the state of the system ceases to change with time.This is a consequence of performing a deterministic analysis, which neglects stochastic effects, and our choice of replication rate in (Equation 10) which decreases with total copy number when w T + δm T > κ and vice versa, guiding the total population to a fixed total copy number.Varying the fission (β) and fusion (γ) rates revealed a negative linear relationship between the steady-state fraction of singletons and copy number (Figure S2B).We may also simulate the system in (Equation 1)-(Equation 9) stochastically, using the stochastic simulation algorithm (Gillespie, 1976), which showed that mean copy number is slightly perturbed from the deterministic prediction due to the influence of variance upon the mean (Grima et al., 2011;Hoitzing, 2017) (Figure 2C).The stationarity of total copy number is a consequence of using δ = 1 for our nominal parametrization (i.e. the line of steady states is also a line of constant copy number).Choosing δ = 1 results in a difference in carrying capacities between the two species, and non-stationarity of mean total copy number, as trajectories spread along the line of steady states to different total copy numbers.Copy number variance initially increases since trajectories are all initialised at the same state, but plateaus because trajectories are constrained in their copy number to remain near the attracting line of steady states (Figure S2C).Mean heteroplasmy remains constant through time under this model (Figure 2D, see (Birky et al., 1983)).This is unsurprising since each species possesses the same replication and degradation rate, so neither species is preferred.
From stochastic simulations we observed that, for sufficiently short times, heteroplasmy variance increases approximately linearly through time for a range of parametrizations (Figure 2E-H), which is in agreement with recent single-cell oocyte measurements in mice (Burgstaller et al., 2018).Previous work has also shown a linear increase in heteroplasmy variance through time for purely genetic models of mtDNA dynamics (see (Johnston and Jones, 2016)).We sought to understand the influence of mitochondrial network dynamics upon the rate of increase of heteroplasmy variance.Wild-type and mutant copy numbers (A) and fused and unfused copy numbers (B) both move towards a line of steady states under a deterministic model, as indicated by arrows.In stochastic simulation, mean copy number (C) is initially slightly perturbed from the deterministic treatment of the system, and then remains constant, while mean heteroplasmy (D) remains invariant with time (see (Equation 75)).In (E)-(H), we show that (Equation 11) holds across many cellular circumstances: lines give analytic results, points are from stochastic simulation.Heteroplasmy variance behaviour is successfully predicted for varying mitophagy rate (E), steady state copy number (F), mutation sensing (G), and fusion rate (H).In (H), fusion and fission rates are redefined as γ → γ0M R and β → β0M where M and R denote the relative magnitude and ratio of the network rates, and γ0, β0 denote the nominal parameterizations of the fusion and fission rates respectively (see Table S1).Figure S2F shows a sweep of M over the same logarithmic range when R = 1.See Figure S2H-P and Table S2 for parameter sweeps numerically demonstrating the generality of the result for different mtDNA control modes.
To this end, we analytically explored the influence of mitochondrial dynamics on mtDNA variability.Assuming that the state of system above (x = (w s , w f , m s , m f )) is initialised at its deterministic steady state (x(t = 0) = x ss ), we took the limit of limit of large mtDNA copy numbers, fast fission-fusion dynamics, and applied a second-order truncation of the Kramers-Moyal expansion (Gardiner, 1985) to the chemical master equation describing the dynamics of the system (see STAR Methods).This yielded a stochastic differential equation for heteroplasmy, via Itô's formula (Jacobs, 2010).Upon forcing the state variables onto the steady-state line (Constable et al., 2016), we derived (Equation 77), which may be approximated for sufficiently short times as x=xss .
(Equation 11) ) is heteroplasmy (a random variable since it is a function of random state variables x), V(h) is the variance of heteroplasmy, µ is the mitophagy rate, n(x) is the total copy number and f s (x) is the fraction of unfused (singleton) mtDNAs, and is thus a measure of the fragmentation of the mitochondrial network.Note that the three functions of state in (Equation 11) have been forced to be their deterministic steady state value (x ss ), which is equivalent to their initial values by model construction.
We find that (Equation 11) closely matches heteroplasmy variance dynamics from stochastic simulation, for sufficiently short times after initialisation, for a variety of parametrizations of the system (Figure 2E-H).
To our knowledge, (Equation 11) reflects the first analytical principle linking mitochondrial dynamics and the cellular population genetics of mtDNA variance.Its simple form allows several intuitive interpretations.As time progresses, replication and degradation of both species occurs, allowing the ratio of species to fluctuate; hence we expect V(h) to increase with time according to random genetic drift (Figure 2E-H).The rate of occurrence of replication/degradation events is set by the mitophagy rate µ, since degradation events are balanced by replication rates to maintain population size; hence, random genetic drift occurs more quickly if there is a larger turnover in the population (Figure 2E).We expect V(h) to increase more slowly in large population sizes, since the birth of e.g. 1 mutant in a large population induces a small change in heteroplasmy (Figure 2F).The factor of h(1 − h) encodes the state-dependence of heteroplasmy variance, exemplified by the observation that if a cell is initialised at h = 0 or h = 1, heteroplasmy must remain at its initial value (since the model above does not consider de novo mutation, see Box 1) and so heteroplasmy variance is zero.Furthermore, the rate of increase of heteroplasmy variance is maximal when a cell's initial value of heteroplasmy is 1/2.In Figure 2G, we show that (Equation 11) is able to recapitulate the rate of heteroplasmy variance increase across different values of δ, which are hypothesized to correspond to different replicative sensing strengths of different mitochondrial mutations (Hoitzing, 2017).We also show in Figures S2D&E that (Equation 11) is robust to the choice of feedback control strength b in (Equation 10).
In (Equation 6), we have made the important assumption that only unfused mitochondria can be degraded via mitophagy, as seen by Twig et al. (2008), hence the total propensity of mtDNA turnover is limited by the number of mtDNAs which are actually susceptible to mitophagy.Strikingly, we find that the dynamics of heteroplasmy variance are independent of the absolute rate of fusion and fission, only depending on the fraction of unfused mtDNAs at any particular point in time (see Figure 2H and Figure S2F).This observation, which contrasts with the model of (Tam et al., 2013(Tam et al., , 2015) ) (see Discussion), arises from the observation that mitochondrial network dynamics are much faster than replication and degradation of mtDNA, by around a factor of β/µ ≈ 10 3 (see Table S1), resulting in the existence of a separation of timescales between network and genetic processes.In the derivation of (Equation 11), we have assumed that fission-fusion rates are infinite, which simplifies V(h) into a form which is independent of the magnitude of the fission-fusion rate.A parameter sweep of the magnitude and ratio of the fission-fusion rates reveals that, if the fusion and fission rates are sufficiently small, (Equation 11) breaks down and V(h) gains dependence upon the magnitude of these rates (see Figure S2H).This regime is, however, for network rates which are approximately 100 times smaller than the biologically-motivated nominal parameterization shown in Figure 2A-D where the fission-fusion rate becomes comparable to the mitophagy rate.Since fission-fusion takes place on a faster timescale than mtDNA turnover, we may neglect this region of parameter space as being implausible.

The influence of mitochondrial dynamics upon heteroplasmy variance under different models of genetic mtDNA control
To demonstrate the generality of this result, we explored several alternative forms of cellular mtDNA control (Johnston and Jones, 2016).We found that when copy number is controlled through the replication rate function (i.e.λ = λ(x), µ = const), when the fusion and fission rates were high and the fixation probability (P (h = 0) or P (h = 1)) was negligible, (Equation 11) accurately described V(h) across all of the replication rates investigated, see Figure S2H-M.The same mathematical argument to show (Equation 11) for the replication rate in (Equation 10) may be applied to these alternative replication rates where a closed-form solution for the deterministic steady state may be written down (see STAR Methods).Interestingly, when copy number is controlled through the degradation rate (i.e.λ =const, µ = µ(x)), heteroplasmy variance loses its dependence upon network state entirely and the f s term is lost from (Equation 11) (see (Equation 86) and Figure S2N-P).A similar mathematical argument was applied to reveal how this dependence is lost (see STAR Methods).
In order to provide an intuitive account for why control in the replication rate, versus control in the degradation rate, determines whether or not heteroplasmy variance has network dependence, we investigated a time-rescaled form of the Moran process (see STAR Methods).The Moran process is structurally much simpler than the model presented above, to the point of being unrealistic, in that the mitochondrial population size is constrained to be constant between consecutive time steps.Despite this, the modified Moran process proved to be insightful.We find that, when copy number is controlled through the replication rate, the absence of death in the fused subpopulation means the timescale of the system (being the time to the next death event) is proportional to f s .In contrast, when copy number is controlled through the degradation rate, the presence of a constant birth rate in the entire population means the timescale of the system (being the time to the next birth event) is independent of f s (see (Equation 98) and surrounding discussion).
A combination of the Moran model and our stochastic results allow us to make a set of biological predictions about interventions to modulate mtDNA behaviour.In Box 1, we use a straightforward Moran model with mutations (infinite sites model) to explore how biological interventions to control network state, mitophagy and copy number influence the rate of accumulation of de novo mutations within a non-dividing cell.Using a Moran model to explore this question is computationally more straightforward than modifying the network model in (Equation 1)-(Equation 10), and more amenable to intuitive understanding, but still effectively captures heteroplasmy variance dynamics (see STAR Methods and Figure S3).We find that reducing both the fraction of singletons and mitophagy rate reduce the rate of accumulation of distinct mutations, as well as the number of mutations per mtDNA.We find that the number of distinct mutations increases with time, as well as the total population size of mtDNAs.We discuss existing evidence in the literature for our model and suggest potential interventions to alleviate the accumulation of mutations through neutral genetic drift.
Box 1. Mathematical modelling of mitochondrial network and genetic processes suggest control strategies against mutant expansions.
In this study, we have argued that the rate of increase of heteroplasmy variance, and therefore the rate of accumulation of pathologically mutated cells within a tissue, increases with mitophagy rate (µ), decreases with total mtDNA copy number per cell (n) and increases with the fraction of unfused mitochondria (termed "singletons", f s ), see (Equation 11).Below, we explore how biological modulation of these variables influences the accumulation of mutations.We use this new insight to propose three classes of strategy to control mutation accumulation and hence address associated issues in aging and disease, and discuss these strategies through the lens of existing biological literature.

TARGETING NETWORK STATE
In (Equation 97) we argue that when population size is controlled in the birth rate, the inter-event rate (Γ) is effectively rescaled by the fraction of unfused mitochondria, i.e.Γ = µnf s .We may use this idea to explore the role of the mitochondrial network in the accumulation of de novo mutations using an infinite sites Moran model (Kimura, 1969) (see Figure B1A).Single cells were modelled over time as having a fixed mitochondrial copy number (n), and at each time step one mtDNA is randomly chosen for duplication and one (which can be the same) for removal.The individual replicated incurs Q de novo mutations, where Q is binomially distributed according to 12) where Binomial(N, p) is a binomial random variable with N trials and probability p of success.L mtDNA = 16569 is the length of mtDNA in base pairs and η = 5.6 × 10 −7 is the mutation rate per base pair per doubling (Zheng et al., 2006); hence each base pair is idealized to have an equal probability of mutation upon replication.
Figure B1B shows that in the infinite sites model, the consequence of (Equation 97) is that the rate of accumulation of mutations per cell reduces as the mitochondrial network becomes more fused, as does the mean number of mutations per mtDNA (Figure B1C).These observations are intuitive: since fusion serves to shield the population from mitophagy, mtDNA turnover slows down, and therefore there are fewer opportunities for replication errors to occur per unit time.Different values of f s in Figures B1B&C therefore correspond to a rescaling of time i.e. stretching of the time-axis.
Concurrently, a study by Chen et al. (2010) observed the effect of deletion of two proteins which are involved in mitochondrial fusion (Mfn1 and Mfn2) in mouse skeletal muscle.Fragmentation of the mitochondrial network induced severe depletion of mtDNA copy number (as found in Figure S2B).The authors also observed that the number of mutations per base pair increased upon fragmentation, which is in agreement with the infinite sites model where fragmentation effectively results in a faster turnover of mtDNA (Figure B1F).
Our models predict that promoting mitochondrial fusion has a two-fold effect: firstly, it slows the increase of heteroplasmy variance (see (Equation 11) and Figure 2H); secondly, it reduces the rate of accumulation of distinct mutations (see Figure B1B&C).These two effects are both a consequence of mitochondrial fusion rescaling the time to the next turnover event, and therefore the rate of random genetic drift.As a consequence, we predict that promoting fusion in youth (assuming mean heteroplasmy is low) could slow down the accumulation and spread of mitochondrial mutations, and perhaps slow aging.Indeed, the premature aging phenotypes of the mitochondrial mutator mouse, which rapidly accumulates mtDNA mutations (Trifunovic et al., 2004), can be rescued through endurance exercise with a resultant normalisation of mitochondrial morphology (Safdar et al., 2011).Since endurance exercise regulates mitochondrial network processes (Yan et al., 2012), it is possible that promotion of mitochondrial fusion contributes to the rejuvenation of the mutator mouse.
If we assume that fusion is selective in favour of wild-type mtDNAs, which appears to be the case at least for some mutations under therapeutic conditions (Suen et al., 2010;Kandul et al., 2016), we predict that a balance between fusion and fission is the most effective means of removing mutant mtDNAs (Figure 3A), perhaps explaining why mitochondrial networks are often observed to exist as balanced between mitochondrial fusion and fission (Sukhorukov et al., 2012;Zamponi et al., 2018).In contrast, if selective mitophagy pathways are induced then promoting fragmentation is predicted to accelerate the clearance of mutants (Figure 3B).

TARGETING MITOPHAGY RATE
Alterations in the mitophagy rate µ have a comparable effect to changes in f s in terms of reducing the rate of heteroplasmy variance (see (Equation 11)) and the rate of de novo mutation (Figure B1B&C) since they both serve to rescale time.Our theory therefore suggests that inhibition of basal mitophagy in youth may be able to slow down the rate of random genetic drift, and perhaps healthy aging, by locking-in low levels of heteroplasmy.Indeed, it has been shown that mouse oocytes (Boudoures et al., 2017) as well as mouse hematopoietic stem cells (de Almeida et al., 2017) have comparatively low levels of mitophagy, which is consistent with the idea that these pluripotent cells attempt to minimise genetic drift by slowing down mtDNA turnover.A previous modelling study has also shown that mutation frequency increases with mitochondrial turnover (Poovathingal et al., 2009).
Alternatively, it has also been shown that the presence of heteroplasmy, in genotypes which are healthy when present at 100%, can induce fitness disadvantages (Acton et al., 2007;Sharpley et al., 2012;Bagwan et al., 2018).In cases where heteroplasmy itself is disadvantageous, especially in later life where such mutations may have already accumulated, accelerating heteroplasmy variance increase to achieve fixation of a species could be advantageous.However, this will not avoid cell-to-cell variability, and the physiological consequences for tissues of such mosaicism is unclear.

TARGETING COPY NUMBER
To investigate the role of mtDNA copy number (mtCN) on the accumulation of de novo mutations, we set f s = 1 such that Γ = µn (i.e. a standard Moran process).We found that varying mtCN did not affect the mean number of mutations per molecule of mtDNA (Figure B1C, inset).However, as the population size becomes larger, the total number of distinct mutations increases accordingly (Figure B1D).In contrast to our predictions, a recent study by Wachsmuth et al. (2016) found a negative correlation between mtCN and the number of distinct mutations in skeletal muscle.However, Wachsmuth et al. (2016) also found a correlation between the number of distinct mutations and age, in agreement with our model.Furthermore, the authors used partial regression to find that age was more explanatory than mtCN in explaining the number of distinct mutations, suggesting age as a confounding variable to the influence of copy number.Our work shows that, in addition to age and mtCN, turnover rate and network state also influence the proliferation of mtDNA mutations.Therefore, one would ideally account for these four variables for jointly, in order to fully constrain our model.
A study of single neurons in the substantia nigra of healthy human individuals found that mtCN increased with age (Dölle et al., 2016).Furthermore, mice engineered to accumulate mtDNA deletions through faulty mtDNA replication (Trifunovic et al., 2004) display compensatory increases in mtCN (Perier et al., 2013), which potentially explains the ability of these animals to resist neurodegeneration.It is possible that the observed increase in mtCN in these two studies is an adaptive response to slow down random genetic drift (see (Equation 11)).In contrast, mtCN reduces with age in skeletal muscle (Wachsmuth et al., 2016), as well as in a number of other tissues such as pancreatic islets (Cree et al., 2008) and peripheral blood cells (Mengel-From et al., 2014).Given the beneficial effects of increased mtCN in neurons, long-term increases in mtCN in youth could delay other age-related pathological phenotypes.

Optimal mitochondrial network configurations for mitochondrial quality control
Whilst the above models of mtDNA dynamics are neutral (i.e.m and w share the same replication and degradation rates), it is often proposed that damaged mitochondria may experience a higher rate of degradation (Narendra et al., 2008;Kim et al., 2007).There are two principal ways in which selection may occur on Figure 3. Selective fusion implies intermediate fusion rates are optimal for mutant clearance whereas selective mitophagy implies complete fission is optimal.Numerical exploration of the shift in mean heteroplasmy for varying fusion/fission ratio, across different selectivity strengths.Stochastic simulations for mean heteroplasmy, evaluated at 1000 days, with an initial condition of h = 0.3 and n = 1000; the state was initialised on the steady state line for the case of f = m = 0, for 10 4 iterations.(A) For selective fusion (see (Equation 13)), for each value of fusion selectivity ( f ), the fusion rate (γ) was varied relative to the nominal parametrization (see Table S1).When f > 0, the largest reduction in mean heteroplasmy occurs at intermediate values of the fusion rate; a deterministic treatment reveals this to be true for all fusion selectivities investigated (see Figure S4).(B) For selective mitophagy (see (Equation 14)), when mitophagy selectivity m > 0, a lower mean heteroplasmy is achieved, the lower the fusion rate (until mean heteroplasmy = 0 is achieved).Hence, complete fission is the optimal strategy for selective mitophagy.
mutant species.Firstly, mutant mitochondria may be excluded preferentially from the mitochondrial network in a background of unbiased mitophagy.If this is the case, mutants would be unprotected from mitophagy for longer periods of time than wild-types, and therefore be at greater hazard of degradation.We can alter the fusion rate (γ) in the mutant analogues of (Equation 1),(Equation 2) and (Equation 7)-(Equation 9) by writing γ → γ/(1 + f ) (Equation 13) for all fusion reactions involving 1 or more mutant mitochondria where f > 0. The second potential selective mechanism we consider is selective mitophagy.In this case, the degradation rate of mutant mitochondria is larger than wild-types, i.e. we modify the mutant degradation reaction to 14) for m > 0.
In these two settings, we explore how varying the fusion rate for a given selectivity ( f and m ) affects the extent of reduction in mean heteroplasmy.Figure 3A shows that, in the context of selective fusion ( f > 0) and non-selective mitophagy ( m = 0) the optimal strategy for clearance of mutants is to have an intermediate fusion/fission ratio.This was observed for all fusion selectivities investigated (see Figure S4) Intuitively, if the mitochondrial network is completely fused then, due to mitophagy only acting upon smaller mitochondrial units, mitophagy cannot occur -so mtDNA turnover ceases and heteroplasmy remains at its initial value.In contrast, if the mitochondrial network completely fissions, there is no mitochondrial network to allow the existence of a quality control mechanism: both mutants and wild-types possess the same probability per unit time of degradation, so mean heteroplasmy does not change.Since both extremes result in no clearance of mutants, the optimal strategy must be to have an intermediate fusion/fission ratio.
In contrast, in Figure 3B, in the context of non-selective fusion ( f = 0) and selective mitophagy ( m > 0), the optimal strategy for clearance of mutants is to completely fission the mitochondrial network.Intuitively, if mitophagy is selective, then the more mtDNAs which exist in fragmented organelles, the greater the number of mtDNAs which are susceptible to selective mitophagy, the greater the total rate of selective mitophagy, the faster the clearance of mutants.

Discussion
In this work, we have sought to unify our understanding of three aspects of mitochondrial physiologythe mitochondrial network state, mitophagy, and copy number -with genetic dynamics.The principal virtue of our modelling approach is its simplified nature, which makes general, analytic, quantitative insights available for the first time.Our models capture the central observation that fused mitochondria are at a lower susceptibility to degradation than unfused mitochondria, since mitophagy acts upon small mitochondrial fragments (Twig et al., 2008).In using parsimonious models, we are able to make the first analytic link between the mitochondrial network state and heteroplasmy dynamics.This is in contrast to other computational studies in the field, whose structural complexity make analytic progress difficult, and accounting for their predicted phenomena correspondingly more challenging.
Using our approach we find that, for a wide class of models, the rate of linear increase of heteroplasmy variance is modulated in proportion to the fraction of unfused mitochondria (see (Equation 11)), which serves to rescale time.This central observation provides a substantial change in our understanding of mitochondrial genetics, as it suggests that the mitochondrial network state, in addition to mitochondrial turnover and copy number, must be accounted for in order to predict the rate of spread of mitochondrial mutations in a cellular population.Crucially, we find that the dynamics of heteroplasmy variance is independent of the absolute rate of fission-fusion events, since network dynamics occur approximately 10 3 times faster than mitochondrial turnover, inducing a separation of timescales.The independence of the absolute rate of network dynamics makes way for the possibility of gaining genetic information about the mitochondrial network experimentally via the network, without the need to quantify absolute fission-fusion rates.By linking with classical statistical genetics, we find that the mitochondrial network also modulates the rate of accumulation of de novo mutations, also due to the fraction of unfused mitochondria serving to rescale time.We find that, in the context of mitochondrial quality control through selective fusion, an intermediate ratio the fusion/fission ratio is optimal due to the finite selectivity of fusion.This latter observation perhaps provides an indication for the reason why we observe mitochondrial networks in an intermediate fusion state under physiological conditions (Sukhorukov et al., 2012;Zamponi et al., 2018).
Our approach has, broadly speaking, been to consider neutral models of mtDNA genetic dynamics.It is, however, typically suggested that increasing the rate of mitophagy promotes mtDNA quality control, and therefore shrinks the distribution of heteroplasmies towards 0% mutant (see (Equation 13) and (Equation 14)).If this is the case, then a neutral genetic model appears to be inappropriate, as mutants experience a higher rate of degradation.Stimulation of the PINK1/Parkin pathway has been shown to select against deleterious mtDNA mutations in vitro (Suen et al., 2010) and in vivo (Kandul et al., 2016), as has repression of the mTOR pathway via treatment with rapamycin (Dai et al., 2013;Kandul et al., 2016).However, the necessity of performing a genetic/pharmacological intervention to clear mutations via this pathway suggests that the ability of tissues to selectively remove mitochondrial mutants under physiological conditions is weak.Consequently, neutral models such as our own are useful in understanding how the distribution of heteroplasmy evolves through time under physiological conditions.Indeed, it has been recently shown that mitophagy is basal (McWilliams et al., 2016) and can proceed independently of PINK1 in vivo (McWilliams et al., 2018), perhaps suggesting that mitophagy has non-selective aspects -although this is yet to be verified conclusively.
We have presented the case of copy number control in the replication rate as being a more intuitive model than control in the degradation rate.The former has the interpretation of biogenesis being varied to maintain a constant population size, with all mtDNAs possessing a characteristic lifetime.The latter has the interpretation of all mtDNA molecules being replicated with a constant probability per unit time, regardless of how large or small the population size is, and changes in mitophagy acting to regulate population size.Such a control strategy seems wasteful in the case of stochastic fluctuations resulting in a population size which is too large, and potentially slow if fluctuations result in a population size which is too small.Furthermore, control in the replication rate means that the mitochondrial network state may act as an additional axis for the cell to control heteroplasmy variance (Figure 2) and the rate of accumulation of de novo mutations (Figure B1B&C).Single-mtDNA tracking through confocal microscopy in conjunction with mild mtDNA depletion could shed light on whether the probability of degradation per unit time per mtDNA varies when mtDNA copy number is perturbed, and therefore provide evidence for or against these two possible control strategies.
Our general approach reveals some apparent differences with previous studies which link mitochondrial genetics with network dynamics.Firstly, Tam et al. (2013Tam et al. ( , 2015) ) found that slower fission-fusion dynamics resulted in larger increases in heteroplasmy variance with time, in contrast to (Equation 11) which only depends on fragmentation state and not absolute network rates.The simulation approach of Tam et al. (2013Tam et al. ( , 2015) ) allowed for mitophagy to act on whole mitochondria, where mitochondria consist of multiple mtDNAs.Faster fission-fusion dynamics tended to form heteroplasmic mitochondria whereas slower dynamics formed homoplasmic mitochondria.It is intuitive that mitophagy of a homoplasmic mitochondrion induces a larger shift in heteroplasmy than mitophagy of a single mtDNA, hence slower network dynamics form more homoplasmic mitochondria.However, this apparent difference with our findings can naturally be resolved if we consider the regions in parameter space where the fission-fusion rate is much larger than the mitophagy rate, as is empirically observed to be the case (Cagalinec et al., 2013;Burgstaller et al., 2014).If the fission-fusion rates are sufficiently large to ensure heteroplasmic mitochondria, then further increasing the fission-fusion rate is unlikely to have an impact on heteroplasmy dynamics.Hence, this finding is potentially compatible with our study, although future experimental studies investigating intra-mitochondrial heteroplasmy would help constrain these models.Tam et al. (2015) also found that fast fission-fusion rates could induce an increase in mean heteroplasmy (Tam et al., 2015), in contrast to Figure 2D which shows that mean heteroplasmy is constant with time after a small initial transient due to stochastic effects.We may speculate that the key difference between our treatment and that of Tam et al. (2013Tam et al. ( , 2015) ) is the inclusion of cellular subcompartments whereby fusion/fission induces migration between such compartments, which may induce dynamics in mean heteroplasmy.We note the uncertainty in accounting for the phenomena observed in such complex models highlights the virtues of a simplified approach which may yield interpretable laws and principles through analytic treatment.
The study of Mouli et al. (2009) suggested that, in the context of selective fusion, higher fusion rates are optimal.This initially seems to contrast with our finding which states that intermediate fusion rates are optimal for the clearance of mutants (Figure 3A).However, the high fusion rates in that study do not correspond directly to the highly fused state in our study.Fission automatically follows fusion in (Mouli et al., 2009), ensuring at least partial fragmentation, and the high fusion rates for which they identify optimal clearing are an order of magnitude lower than the top fusion rate they consider.In the case of complete fusion, mitophagy cannot occur in the model of Mouli et al. (2009), so there is no mechanism to remove dysfunctional mitochondria.It is perhaps more accurate to interpret the observations of Mouli et al. (2009) as implying that selective fusion shifts the optimal fusion rate higher, when compared to the case of selective mitophagy alone.Therefore the study of Mouli et al. (2009) is compatible with Figure 3A.Furthermore, Mouli et al. (2009) also found that when fusion is non-selective and mitophagy is selective, intermediate fusion rates are optimal whereas Figure 3B shows that complete fragmentation is optimal for clearance of mutants.Optimality of intermediate fusion in the context of selective mitophagy in the model of Mouli et al. (2009) likely stems from two aspects of their model: i) mitochondria consist of several units which may or may not be functional; ii) the sigmoidal relationship between number of functional units per mitochondrion and mitochondrial 'activity' (the metric by which optimality is measured).Points (i) and (ii) imply that small numbers of dysfunctional mitochondrial units have very little impact on mitochondrial activity, so fusion may boost total mitochondrial activity in the context of small amounts of mutation.So whilst Figure 3B remains plausible in light of the study of (Mouli et al., 2009) if reduction of mean heteroplasmy is the objective of the cell, it is also plausible that non-linearities in mitochondrial output under cellular fusion (Hoitzing et al., 2015) result in intermediate fusion being optimal in terms of energy output in the context of non-selective fusion and selective mitophagy.Future experimental studies quantifying the importance of selective mitophagy under physiological conditions would be beneficial for understanding heteroplasmy variance dynamics.We note again that the ubiquity of heteroplasmy (Payne et al., 2012;Ye et al., 2014;Morris et al., 2017) suggests that a neutral drift approach to mitochondrial genetics may be justified, which contrasts with the studies of Tam et al. (2013Tam et al. ( , 2015) ) and Mouli et al. (2009) which focus purely on the selective effects of mitochondrial networks.
In order to fully test our model, further single-cell longitudinal studies are required.For instance, the study by Burgstaller et al. (2018) found a linear increase in heteroplasmy variance through time in single oocytes.
Our work here has shown that measurement of the network state, as well as turnover and copy number, are required to account for the rate of increase in heteroplasmy variance.Joint longitudinal measurements of f s , µ and n, with heteroplasmy quantification, would allow verification of (Equation 11) and aid in determining the extent to which neutral genetic models are explanatory.This could be achieved, for instance, using the mito-QC mouse (McWilliams et al., 2016) which allows visualisation of mitophagy and mitochondrial architecture in vivo.Measurement of f s , µ and n, followed by e.g.destructive single-cell whole-genome sequencing of mtDNA would allow validation of how µ, n and f s influence V(h) and the rate of de novo mutation (see Figure B1).One difficulty is sequencing errors induced through e.g.PCR, which hampers our ability accurately measure mtDNA mutation within highly heterogeneous samples (Woods et al., 2018).Morris et al. (2017) have suggested that single cells are highly heterogeneous in mtDNA mutation, with each mitochondrion possessing 3.9 single-nucleotide variants on average.Error correction strategies during sequencing may pave the way towards high-accuracy mtDNA sequencing (Woods et al., 2018;Salk et al., 2018), and allow us to better constrain models of heteroplasmy dynamics.

CONTACT FOR REAGENT AND RESOURCE SHARING
As Lead Contact, Iain Johnston is responsible for all requests for further information.Please contact, Iain Johnston (i.johnston.1@bham.ac.uk) with requests and inquiries.

METHOD DETAILS
Constant rates yield unstable copy numbers for a model describing mtDNA genetic and network dynamics We explored a simpler network system than the one presented in the Main Text, but found that it produced instability in mtDNA copy numbers, which we regard as biologically undesirable.Consider the following set of Poisson processes for singleton (s) and fused (f ) species 15) 16) 17) 18) 19) 20) 21) where (Equation 15)-(Equation 17) are analogous to (Equation 1)-(Equation 3) where mutant species are neglected.The parameter ρ is shared amongst all of the birth and death reactions in (Equation 18)-(Equation 21).ρ represents the intuitive assumption that, in order for a stable population size to exist, birth should balance death.However, for the network to have any effect at all, singletons should be at an increased risk of mitophagy relative to fused species.We represent the increased risk of singleton mitophagy with the parameter η.Since additional death is introduced into the system when η > 1, we include the parameter α > 1 as an increased global biogenesis rate to balance the increased mitophagy of singletons.We may write the above system as a set of ordinary differential equations ds dt = −γs 2 − γf s + βf + αρs − ηρs (Equation 22) 23) where we have enforced the stochastic reaction rate to be equivalent to the deterministic reaction rate, and hence the s 2 term is proportional to γ rather than 2γ (justification of this is presented below, see (Equation 34)).In Figure S1 we see that the system displays a trivial steady state at s = f = 0 and a non-trivial steady state.Computing the eigenvalues of the Jacobian matrix at the non-trivial steady state indicates that it is a saddle node, and therefore unstable.Initial copy numbers which are too small tend towards extinction with time, and initial copy numbers which are too large tend towards a copy number explosion.This simple example suggests that a system of this form with constant reaction rates is unstable, and therefore biologically unlikely to exist under reasonable circumstances.We hence consider analogous biochemical reaction networks with a replication rate which is a function of state, to prevent extinction and divergence of the total population size.

Conversion of a chemical reaction network into ordinary differential equations
The following section outlines the steps in converting a set of chemical reactions into a set of ordinary differential equations (ODEs).In particular, we pay special attention to the fact that the rate of a chemical reaction with a stochastic treatment is not always equivalent to the rate in a deterministic treatment (Wilkinson, 2011), as we will explain below.This subtlety is sometimes overlooked in the literature.This section draws on a number of standard texts (Gillespie, 1976;Van Kampen, 1992;Gillespie, 2007;Wilkinson, 2011) as well Grima (2010).We hope this harmonized treatment will be of help as a future reference.
Consider a general chemical system consisting of N distinct chemical species (X i ) interacting via R chemical reactions, where the j th reaction is of the form 24) where s ij and r ij are stoichiometric coefficients.We define kj as the microscopic rate for this reaction.The dimensionality of this parameter will vary depending upon the stoichiometric coefficients s zj .kj may be loosely interpreted as setting the characteristic timescale (i.e. the cross section (Wilkinson, 2011)) of reaction j.
The chemical master equation (CME) describes the dynamics of the joint distribution of the state of the system and time, moving forwards through time.Defining the state of the system as x = (x 1 , . . ., x N ) T , where x i is the copy number of the i th species, allows us to write the CME as (Grima, 2010) 25) where Ω is the volume of the compartment in which the reactions occur (also known as the system size), S ij = r ij − s ij is the stoichiometry matrix, and E is referred to as the step operator and is defined through the relation E −Sij i (g(x)) = g(x 1 , . . ., x i − S ij , . . ., x N ), for any function of state g(x).fj (x, Ω) is the microscopic rate function of reaction j, which in general depends on both the state and the system size.A factor of Ω is explicitly included in this definition of the chemical master equation so that our treatment is compatible with Van Kampen's system size expansion (Van Kampen, 1992).As a consequence of this, the probability that, given the current state x, the j th reaction occurs in the time interval [t, t + dt) somewhere in Ω (Gillespie, 2007) is âj (x, Ω)dt := Ω fj (x, Ω)dt.
For the microscopic rate function, we may write fj (x, Ω) = kj 27) This equation counts the number of available combinations of reacting molecules (Gillespie, 1976;Wilkinson, 2011), whilst taking into account scaling with system size (Grima, 2010).
We also introduce the deterministic rate equation (generally considered to be the macroscopic analogue of the CME) which is defined as (Van Kampen, 1992;Grima, 2010) 28) where φ = (φ 1 , . . ., φ N ) T is the vector of macroscopic concentrations (of dimensions molecules per unit volume) and fj (φ) is the macroscopic rate function satisfying 29) where kj is the macroscopic rate for the j th reaction.Note that we distinguish between kj and kj , respectively the rate constants for the discrete and continuous pictures, although this distinction is sometimes not emphasized in the literature (Grima, 2010;Grima et al., 2011;Van Kampen, 1992).The physical meaning of kj is not immediately obvious: we argue that this parameter only gains physical meaning through the following procedure.
As stated by Wilkinson (2011), if we intend for the microscopic description in (Equation 25) to correspond to the macroscopic description in (Equation 28), the rate of consumption/production of particles for every reaction must be the same in the deterministic limit of the stochastic system (the conditions for which we define below).Therefore, we apply the following constraint in the limit of large copy numbers lim xi→∞ fj (x, Ω) = fj (φ) ∀ i, j.
(Equation 30) In applying this constraint on all species i and all reactions j, we may derive a general relationship between kj and kj (Equation 31) We can make two approximations to generate a more convenient relationship between the microscopic and macroscopic rates.Firstly, we assume that 32) This is a small noise approximation, since it is often assumed that x i = Ωφ i + Ω 1/2 ξ i , where ξ i is a noise term (Van Kampen, 1992).If ξ i is small then x i ≈ Ωφ i is a valid approximation.Secondly, we assume that 33) This is a large copy number approximation: in the case of e.g. a bimolecular reaction (2X i → * ) with s ij = 2, the approximation is of the form x i (x i − 1) ≈ x 2 i or x i ≈ x i − 1.By applying (Equation 33) to the factor of x i ! in (Equation 31), the factor of (x i − s ij )! cancels from the left-hand side.Simplifying using (Equation 32), φ sij i cancels from both sides and we arrive at the important relationship 34) With (Equation 28), (Equation 29) and (Equation 34) one may therefore write down a set of ODEs for an arbitrary chemical reaction network, with constant reaction rates, in terms of the microscopic rates kj .This equation highlights that for reactions with s ij ≥ 2, kj = kj , as is the case for bimolecular reactions of the form 2X i → * (see (Equation 1) and (Equation 15)).
It is important to note that if the microscopic rate function is a function of state then k = k(x) and k = k(φ) ≈ k(x/Ω).In this case, (Equation 34) still applies since the above argument assumed nothing about the particular forms of k and k.However, additional factors of Ω −1 are induced by applying (Equation 32), which may carry through to the individual parameters of k(φ).A demonstration of this is given in the following section.

Deriving an ODE description of the mitochondrial network system
In this section we show how to derive an ODE description of the network system described in (Equation 1)-(Equation 9) in the Main Text.In accordance with the notation in the previous section, we will redefine all of the rates in (Equation 1)-(Equation 10) with a hat notation (â, for a general rate parameter a), to reflect that these are stochastic rates.Deterministic rates will be denoted with a tilde (ã).Our aim will be to write a set of ODEs in terms of the stochastic rates, â, for which we are able to estimate values.
We will begin by considering the fusion network equations (Equation 1) and (Equation 2).For clarity, we rewrite (Equation 1) to allow the reaction to proceed with some arbitrary rate ρ: 35) where X denotes either a wild-type (W ) or mutant (M ).We will subsequently fix ρ to the rate of all other fusion reactions γ.We do this because (Equation 1) is a bimolecular reaction involving one species: a fundamentally different reaction to bimolecular reactions involving two species, as we will now see.Since ρ, γ = const, we may use (Equation 34), resulting in the deterministic rates ρ = ρ 2 (Equation 36) γ = γ (Equation 37) for (Equation 35) and (Equation 2) respectively.If we then enforce the microscopic rates to be equal for both of these fusion reactions, i.e. ρ = γ, then ρ = γ/2.All other fusion reactions have γ = γ by application of (Equation 34).Application of (Equation 34) to the fission reaction in (Equation 3) shows that β = β.For (Equation 4), we have chosen a λ which is not a constant, but a function of the copy numbers of the chemical species ( λ = λ(x) where x = (w s , w f , m s , m f ), see (Equation 10)).As pointed out in the previous section, care must be taken in writing down the deterministic analogue of λ.Applying (Equation 34), we have (Equation 38) Applying (Equation 32) and equating individual terms, we arrive at μ = μ (Equation 39) bw s = b w s Ω =⇒ b = b/Ω (Equation 40) bk = bk =⇒ k = kΩ (Equation 41) bδ m s = bδ φ ms =⇒ δ = δ.
(Equation 42) In this study, we let Ω = 1 so the above 4 parameters are identical to their deterministic counterparts.Hence, by application of (Equation 28), we arrive at the following set of ODEs (Equation 43) 44) 45) 46) The steady state solution of this system of ODEs may be calculated, but its form is complex.For notational simplicity, we will drop the hat notation.Defining (Equation 47) the non-trivial, physically-realizable, component of the steady state may be parameterised in terms of φ ms and written as (Equation 48) 49) (Equation 50) Note that, since the steady state is parametrized by φ ms , the steady state is therefore a line.

Proof of heteroplasmy relation for linear feedback control
In this section we show that (Equation 11) holds for the system described by (Equation 1)-(Equation 9) given the replication rate in (Equation 10) using the Kramers-Moyal expansion under conditions of large copy number and fast network churn (to be defined below); the approach used here is similar to Constable et al. (2016).Consonant with the self-contained objectives of STAR methods, we draw together elements from the literature to provide a coherent derivation; we therefore hope that the following exposition may provide clarity for a wider audience.
Kramers-Moyal expansion of the chemical master equation for large copy numbers Customarily, the Kramers-Moyal expansion is formed using a continuous-space notation (Gardiner, 1985), so we will initially proceed in this way.Following the treatment by Gardiner (1985), we begin by re-writing the chemical master equation (Equation 25) (CME) as 51) where we have set Ω = 1.T (x|x ) is the transition rate from state x → x, and the dependence upon the initial condition has been suppressed for notational convenience.We now proceed by expanding the CME.
The multivariate Kramers-Moyal expansion may be written as 52) where H(x) is the Hessian matrix of T (x |x)P (x) 53) (see (Gardiner, 1985) for a proof of this in the univariate case).
Note that a transition to each possible neighbouring state x corresponds to some reaction j which moves the state from x → x .Since we know the influence of each reaction on state x through the constant stoichiometry matrix S ij , and that the propensity of a reaction does not depend upon x itself (see (Equation 27)), we may transition from a notation involving x and x into a notation involving x and j.We may therefore define T j (x) := T (x |x) ≡ fj (x) (see (Equation 27)), and let H(x) → H j (x).
We now make a large copy number assumption in order to simplify T j (x).To take a large copy number limit, we assume that 54) Note that this approximation is exact when s ij = 0, 1, but inexact when s ij ≥ 2. For example, if we consider the second-order bimolecular reaction in (Equation 1), (Equation 54) is equivalent to assuming w 2 s ≈ w s (w s − 1); consequently, a factor of 1/(s ij ! ) = 1/2 arises in T j (x) as a combinatorial factor from stochastic considerations.

Fokker-Planck equation for chemical reaction networks
We now wish to re-write (Equation 52) as a Fokker-Planck equation.Since the integral in (Equation 52) is over x , and every x corresponds to a reaction j, we may interpret the integral in (Equation 52) as a sum over all reactions, i.e. dx → R j=1 .Hence, for the j th reaction, [(x − x)] i = S ij .With these observations, we may write the first integral of (Equation 52) as 55) where 56) 24), and T is the vector of transition rates, of length R (for which we have taken a large copy number approximation in (Equation 54)).To re-write the second integral of (Equation 52), we write an element of the Hessian H j in (Equation 53) as 57) where j = 1, . . ., R and l, m = 1, . . ., N .Thus, we may write 58) where B is an N × N matrix, and Diag(Y) is a diagonal matrix whose main diagonal is the vector Y.We may therefore re-write (Equation 52) as a Fokker-Planck equation for the state vector x of the form (Equation 60) Fokker-Planck equation for an arbitrary function of state We now wish to make a change of variables in (Equation 60) to write down a Fokker-Planck equation for an arbitrary scalar function of state x (which we will later set to be heteroplasmy).To do this, we wish to make use of Itô's formula, which allows a change of variables for an SDE.In general, the Fokker-Planck equation in (Equation 60) is equivalent (Jacobs, 2010) to the following Itô stochastic differential equation (SDE) dx = A dt + G dW (Equation 61) where GG T ≡ B and dW is a vector of independent Wiener increments of length N , and a Wiener increment dW satisfies 2t) .(Equation 62) Itô's formula states that, for an arbitrary function h(x, t) where x satisfies (Equation 61), we may write the following SDE dh(x, t) = ∂h ∂t + (∇h) where H h (x) is the Hessian matrix of h(x, t) (see (Equation 53), where T (x |x)P (x) should be replaced with h(x, t)).Given the form of B in (Equation 59) we let which satisfies GG T ≡ B.
For convenience, we may also perform the transformation purely at the level of Fokker-Planck equations.Let h(x, t) satisfy the general Fokker-Planck equation 65) for scalar functions Ã(h, t) and B(h, t).Using the cyclic property of the trace in (Equation 63), we may identify Ã = ∂h ∂t + (∇h) 66) where Tr is the trace operator.Also, from (Equation 63), (Equation 67) Hence, using (Equation 65), (Equation 66) and (Equation 67), we may write down a Fokker-Planck equation for an arbitrary function of state in terms of A and B.
An SDE for heteroplasmy forced onto the steady state line in the high-churn limit It has been demonstrated that SDE descriptions of stochastic systems which possess a globally-attracting line of steady states may be formed in the long-time limit by forcing the state variables onto the steady state line (Constable et al., 2016;Parsons and Rogers, 2017).Such descriptions may be formed in terms of a parameter which traces out the position on the steady state line, hence reducing a high-dimensional problem into a single dimension (Constable et al., 2016;Parsons and Rogers, 2017).In our case, heteroplasmy is a suitable parameter to trace out the position on the steady state line.We seek to use similar reasoning to verify (Equation 11).In what follows, we will assume that x(t = 0) = x ss , where x ss is the state which is the solution of A = 0 (which is equivalent to finding the steady state solution of the deterministic rate equation in (Equation 28) due to our assumption of large copy numbers and Ω = 1), so that we may neglect any deterministic transient dynamics.
Inspection of the steady state of the ODE description of our system reveals that the set of steady state solutions forms a line (see (Equation 48)-(Equation 50)).Inspection of the steady state solution reveals that the steady state depends on the fusion (γ) and fission (β) rates.Mitochondrial network dynamics occur on a much faster timescale than the replication and degradation of mtDNA: the former occurring on the timescale of minutes (Twig et al., 2008) whereas the latter is hours or days (Johnston and Jones, 2016).We seek to use this separation of timescales to arrive at a simple form for V(h).We redefine the fusion and fission rates such that 68) where M is a constant which determines the magnitude of the fusion and fission rates, which we call the "network churn".We now wish to use heteroplasmy as our choice for the function of state in the Fokker-Planck equation in (Equation 65).We will first compute the diffusion term B for heteroplasmy using (Equation 67).If we constrain the state x to be forced onto the steady state line x ss (as per (Constable et al., 2016;Parsons and Rogers, 2017)) in the high-churn limit, then upon defining 70) we have lim (Equation 71) (Equation 71) is difficult to understand.In order to perform further simplification, we make an ansatz for the form of B ( BAn ) and seek to determine whether our ansatz is equivalent to the derived form of B under the constraints defined on the left-hand side of (Equation 71).Our ansatz takes the form x=xss (Equation 72) where f s (x) := (w s + m s )/(w s + w f + m s + m f ) and n(x) := w s + w f + m s + m f .Notice that this ansatz is more general than (Equation 71), since it has no explicit dependence upon the parameters of the control law assumed in (Equation 10), and only explicitly depends upon functions of state x.
Upon substituting the steady state x ss into the ansatz in (Equation 72) and taking the high-churn limit, we find that (Equation 73) After some algebra (see GitHub repository for Mathematica notebook), it can be shown that (Equation 71) and (Equation 73) are equivalent, i.e.
BAn ≡ lim (Equation 74) As such, we may use B and BAn interchangeably in the limit of high network churn.Furthermore, it can be shown after some algebra that the drift of heteroplasmy when forced onto the steady state line is 0, i.e.
Ã x=xss ≡ 0. (Equation 75) A similar result is shown in (Constable et al., 2016) ((Equation S59) therein).Substituting h(x, t) = h, Ã and B into the Fokker-Planck equation for an arbitrary function of state (Equation 65), we have 76) which is equivalent to the following SDE for heteroplasmy dW (Equation 77) in the limit of large network churn, large copy numbers, and a second-order truncation of the Kramers-Moyal expansion.Note that, although the state has been forced onto the steady state, stochastic fluctuations mean that trajectories may move along the line of steady states, so the diffusion coefficient is not constant in general.
We may calculate the new value of x ss (h) for every displacement due to Wiener noise in h, and substitute into f s (x) and n(x) to determine the diffusion coefficient at the next time step.However, for sufficiently short times, and large copy numbers (i.e.low diffusivity of h), we may assume that the diffusion coefficient in (Equation 77) may be approximated as constant.Since the general solution of the SDE dy = √ B dW (Equation 78) for B = const is y ∼ N (y|y 0 , Bt) (Equation 79) where N (y|y 0 , σ 2 ) is a Gaussian distribution on y with mean y 0 and variance σ 2 , and y 0 = y(t = 0).Since we have assumed that the state is initialised at x(t = 0) = x ss , there are no deterministic transient dynamics, so we may write 80) where V returns variance of a random variable.Note that in this equation, we take x = x ss = const, since we have assumed a low-diffusion limit.Also note that this equation is of precisely the same form as (Equation 12) of Johnston and Jones (2016), except with an additional proportionality factor of f s induced by the inclusion of a mitochondrial network.

Heteroplasmy variance relations for alternative model structures and modes of genetic mtDNA control
Here we explore the implications of alternative model structures upon (Equation 77).Firstly, we may consider replacing (Equation 4) with (Equation 81) This corresponds to the case where replication coincides with fission, see (Lewis et al., 2016).Repeating the calculation in the previous section also results in (Equation 77), so the result is robust to the particular choice of mtDNA replication reaction (see GitHub repository for Mathematica notebook).Secondly, we may explore the impact of allowing non-zero mtDNA degradation of fused species.This could correspond to autophagy-independent degradation of mtDNA, for example via the exonuclease activity of POLG (Medeiros et al., 2018).To encode this, we may add the following additional reaction 82) where 0 ≤ ξ ≤ 1.We were not able to make analogous analytical progress in this instance.However, numerical investigation (Figure S2G) revealed that the following ansatz was able to predict heteroplasmy variance dynamics 83) In other words, allowing degradation of fused species results in a linear correction to our heteroplasmy variance formula in (Equation 11).If fused species are susceptible to degradation at the same rate as unfused species (ξ = 1), then V(h) loses f s dependence entirely and the mitochondrial network has no influence over heteroplasmy dynamics.We also explored various different forms of λ(x) and µ(x), which we label A-G after (Johnston and Jones, 2016), and X-Z for several newly-considered functional forms, see Table S2 and Figure S2H-P.Note that control D of (Johnston and Jones, 2016) involves no feedback, which we do not explore -see Figure S1, and the discussion surrounding (Equation 15).The argument presented in the previous section requires the steady state solution of the system to be solvable, since we require the explicit form of x ss in (Equation 71), (Equation 73) and (Equation 75).For controls B, C, E, F, G, Y and Z in Table S2, the steady states are solvable and similar arguments to the above can be applied (see the GitHub repository for Mathematica notebooks).Controls B, C, E, F all satisfy (Equation 77); this can be shown numerically for controls A and X.However, controls G, Y and Z satisfy dh dW.
(Equation 84) Notably, (Equation 84) does not depend on f s , unlike (Equation 77) (see GitHub repository for Mathematica notebooks).This is because control of copy number occurs in the degradation rate, rather than the replication rate, for controls G, Y and Z.A modified version of a Moran process (presented below) can provide intuition for why the diffusion rate of heteroplasmy variance depends on the network state when the population is controlled through replication, and does not depend on network state when the population is controlled through degradation.

Choice of nominal parametrization
In this section we discuss our choice of nominal parametrization for the network system in (Equation 1)-(Equation 9), given the replication rate in (Equation 10).We will first discuss our choice of network parameters.Cagalinec et al. (2013) found that the average fission rate in cortical neurons is 0.023±0.003fissions/mitochondria/min. Assuming that this value is representative of the fission rate in general, and converting this to units of per day, we may write the mitochondrial fission rate as β = 33.12day −1 .
Note that the dimensions of β are day −1 and not mitochondrion −1 day −1 .This is because if the propensity (see (Equation 26), where Ω = 1) of e.g.(Equation 3) is âfis,w = βw f then the mean time to the next event is 1/(βw f ); therefore the dimension of β is per unit time and copy numbers are pure numbers, i.e. dimensionless.Similar reasoning constrains the dimension of the fusion rate, see below.
Evaluation of the fusion rate is more involved, since fusion involves two different chemical species coming together to react whereas fission may be considered as spontaneous.Furthermore, there are 7 different fusion reactions whereas there are only 2 fission reactions.For simplicity, assume that all species have a steady-state copy number of x i = 250 (resulting in a total copy number of 1000, heteroplasmy of 0.5 and 50% of mitochondria existing in the fused state).Neglecting subtleties relating to bimolecular reactions involving one species (see (Equation 34)), each fusion reaction proceeds at rate âfus,j ≈ γx 2 i .Since there are 7 fusion reactions ((Equation 1), (Equation 2), (Equation 7)-(Equation 9)), the total fusion propensity is âfus ≈ 7γx 2 i .Similarly, the total fission propensity is âfis = β(w f + m f ) = 2βx i .Since we expect macroscopic proportions of both fused and fissioned species in many physiological settings, we may equate the fusion and fission propensities, âfus = âfis , and rearrange for the fusion rate γ to yield γ = 2βx i /(7x 2 i ) ≈ 3.8 × 10 −2 day −1 .Note that the orders of magnitude difference between β and γ stems from the observation that fusion propensity depends on the square of copy number whereas the fission propensity depends on copy number linearly.
Given the network parameters, we then explored appropriate parametrizations for the genetic parameters: the mitophagy rate (µ) and the parameters of the linear feedback control (κ, b and δ, see (Equation 10)).mtDNA half-life is observed to be highly variable: in mice this can be between 10-100 days (Burgstaller et al., 2014).For consistency with another recent study investigating the relationship between network dynamics and heteroplasmy, we use an mtDNA half-life of 30 days (Tam et al., 2015).
The parameter δ in the replication feedback control (see (Equation 10)) may be interpreted as the "strength of sensing of mutant mtDNA" in the feedback control (Hoitzing et al., 2017).Assuming that fluctuations in copy number of mutants and wild-type molecules are sensed identically (as may be the case for e.g.non-coding mtDNA mutations) we may reasonably assume a model of δ = 1 as the simplest case of a neutral mutation.
We are finally left with setting the parameters κ and b in the linear feedback control (Equation 10).In the absence of a network state and mutants, κ is precisely equal to the steady state copy number, since the degradation rate equals the replication rate when κ = w.However, the presence of a network means that a subpopulation of mtDNAs (namely the fused species) are immune to death, resulting in κ no longer being equivalent to the steady state copy number.The parameter b may be interpreted as the feedback control strength, which determines the extent to which the replication rate changes given a unit change in copy number.
Given a particular value of b, we may search for a κ which gives a total steady state copy number (n) which is closest to some target value (e.g.1000 as a typical total mtDNA copy number per cell in human fibroblasts (Kukat et al., 2011)).We swept a range of different values of b and found that, for values of b smaller than a critical value (b * ), a κ could not be found whose deterministic steady state was sufficiently close to n = 1000.This result is intuitive because in the limit of b → 0, λ = const.From the analysis above we have shown that constant genetic rates (µ, λ) result in unstable copy numbers, and therefore a sufficiently small value of b is not expected to yield a stable non-trivial steady state solution.We chose b ≈ b * , and the corresponding κ, such that the steady state copy number is controlled as weakly as possible given the model structure.

Rate renormalization
In (Equation 1)-(Equation 10) we have neglected reactions such as 85) because they do not change the number of molecules in our state vector x = (w s , w f , m s , m f ).One may ask whether neglecting such reactions means that it is necessary to renormalize the fission-fusion rates which were estimated in the preceding section.In estimating the nominal parametrization above, we began by using a literature value for the mitochondrial fission rate, and then matched the fusion rate such that the summed hazard of a fusion event approximately balanced the fission rate.This matching procedure is reasonable, since we observe a mixture of fused and fissioned mitochondria under physiological conditions: choosing a fusion rate which is vastly different results in either a hyperfused or fragmented network.We must therefore only justify the fission rate.(Equation 3) assumes that a fission reaction always results in a singleton, and a singleton is by definition a molecule which is susceptible to mitophagy (see (Equation 6)).Therefore, if fission reactions always result in mitochondria containing single mtDNAs which are susceptible to mitophagy, then we expect our model to match well to true physiological rates.If, on the other hand, fission reactions between large components of the network which are too large to be degraded are common, then renormalization of β by the fraction of fission events which result in a sufficiently small mitochondrion would be necessary, which would in turn renormalize γ through our matching procedure.We are not aware of experimental measurements of the fraction of fission events which result in mitochondria which contain a particular number of mtDNAs.Such an experiment, combined with the distribution of mitochondrial sizes which are susceptible to mitophagy, would allow us to validate our approach.Despite this, the robustness of our results over approximately 4 orders of magnitude for the fission-fusion rate (Figure S2H-P) provides some indication that our results are likely to hold in physiological regions of parameter space.
A modified Moran process may account for the alternative forms of heteroplasmy variance dynamics under different models of genetic mtDNA control We sought to gain insight into why control of population size through the replication rate (λ = λ(x), µ = const) results in heteroplasmy variance depending on the fraction of unfused mitochondria (see (Equation 11)), whereas control of population size through the degradation rate (µ = µ(x), λ = const) results in heteroplasmy variance becoming independent of network state, where 86) We will proceed by considering an analogous Moran process to the set of reactions presented in (Equation 1)-(Equation 9).First, consider a haploid biallelic Moran process consisting of wild-types and mutants, in a population of fixed size n.At each step in discrete time, a member of the population is chosen for birth, and another for death.Let m t denote the copies of mutants at time t.Then, (Equation 87) 88) Defining h t := m t /n then from (Equation 87) 89) and therefore 90) Suppose that, instead of the process occurring with discrete time, instead the process occurs with continuous time, where each event is a simultaneous birth and death, and is modelled as a Poisson process.Suppose that the next death event) is proportional to f s .We note that this argument is only a heuristic, since the Moran process is defined such that birth and death events occur simultaneously and therefore do not possess separate propensities (Γ birth and Γ death ).
The second case we consider is when each individual has a constant rate of birth, hence Γ birth ∝ n.Then the death rate is chosen such that Γ birth = Γ death .In this case Γ = λn, where λ is a proportionality constant.The same argument from (Equation 91) to (Equation 93) may be applied, with an appropriate rescaling of time, and we arrive at (Equation 86).This is analogous to when µ = µ(x) and λ = const in the network system.Hence, when µ = µ(x) and λ = const, the presence of a constant birth rate in the entire population means the timescale of the system (being the time to the next birth event) is independent of f s .

QUANTIFICATION AND STATISTICAL ANALYSIS
In Figures S2D,H-P, we compare (Equation 11) and (Equation 86) to stochastic simulations, for various parameterizations and replication/degradation rates.In order to quantify the accuracy of these equations in predicting V(h, t), we define the following error metric = 1 − V(h, t) Th E t ( V(h, t) Sim ) (Equation 99) where V(h, t) is the time derivative of heteroplasmy variance with subscripts denoting theory (Th) and simulation (Sim).An expectation over time (E t ) is taken for the stochastic simulations, whereas V(h, t) is a scalar quantity for (Equation 11) and (Equation 86).15)-(Equation 21).The system displays two steady states: a trivial steady state at s + f = 0, and a non-trivial steady state.At a steady state, both time derivatives in (Equation 22) and (Equation 23) vanish.Trajectories (blue) show the evolution of the system up to t = 1000.The direction and magnitude of the derivative at points in space are shown by red arrows.Trajectories can be seen either decaying to s = f = 0 or tending to infinity.γ = 0.01656, β = 33.12,ρ = 0.023, η = 1.1, α = 1.04.
Table S1.Related to Figure 2 and STAR Methods.Nominal parametrization for network system.Dimensions of parameters are derived using individual terms in (Equation 43)-(Equation 46); copy numbers of particular species are pure numbers and therefore dimensionless.Since we have chosen a system size of Ω = 1 throughout, we set the dimension of volume to 1. See STAR Methods for further details of parameter justification.S2 and Method details. (A) Deterministic dynamics of total copy number, which is controlled to a particular steady state value.(B) Defining a knock-down (KD) factor (k −1 = 0.1, 0.2 . . ., 1.0), the fission rate was rescaled to β → β/k (red) and the fusion rate to γ → γ/k (blue), causing a linear increase and decrease in total copy number respectively under a deterministic treatment.(C) Copy number variance for stochastic simulations initially increases, since all stochastic simulations begin with the same initial condition, but then plateaus since the steady state line is globally attracting.(D) Error in (Equation 11) in a sweep over the feedback control strength, b.Dotted line denotes a 5% error according to (Equation 99).(E) V(h) profile for the parameterization with the largest error in (D).(F) Sweeps of the network rate magnitude (see Figure 2H).Heteroplasmy variance is approximately independent of absolute network rates over a broad range of network magnitudes.(G) Allowing fused species to be degraded with relative rate ξ (Equation 82), stochastic simulations for heteroplasmy variance (markers) and (Equation 83) (lines) are shown.Fused species degradation induces a linear correction to the heteroplasmy variance formula.(H-P) Error in (Equation 11) (ansatz) and (Equation 86) (ansatz indep network) for sweeps of the fusion and fission rates the corresponding feedback control functions in Table S2.Equations are accurate to at least 5% (blue regions) across large regions of parameter space, for many control laws.Fusion and fission rates are redefined as γ → γ0M R and β → β0M where M and R denote the magnitude and ratio of the network rates, and γ0, β0 denote the nominal parameterizations of the fusion and fission rates respectively.Summary statistics for 10 4 realizations, with initial condition h = 0.3 and evaluated at t = 500 days.Errors in V(h) (see (Equation 99)) smaller than 5% are truncated and are shown as blue.Parameterizations where a deterministic steady state could not be found for an initial condition of h = 0.3 are shown in grey.Inset figures, where present, display the probability of fixation at h = 0.Where insets are not present, the probability of fixation is negligible.(H-M) When λ = λ(x) and P (h = 0) is low, (Equation 11) performs well in the high-churn limit.(N-P) When µ = µ(x) and P (h = 0) is low, (Equation 86) performs well in the high-churn limit.

Parameter Description
Table S2.Relating to Figure 2 and Figure S2.Nominal parameterizations for the alternative feedback control functions explored.In all cases, the nominal fission and fusion rates were β = 33.12,γ = 0.038 respectively.S1 µ; see Table S1 Control E in (Johnston and Jones, 2016) and GitHub repository

Interpretation
Relaxed replication (see Figure S2I)   13), according to an ODE An ODE treatment allows smaller heteroplasmy changes to be probed without the need to resort to an infeasible number of stochastic simulations.Displaying the relative change in heteroplasmy (∆h) after t = 1000 days.We observe that a reduction in heteroplasmy is achieved at intermediate fusion rates at all non-zero fusion selectivities investigated.Grey denotes a change which is smaller than floating point precision.

Figure 2 .
Figure 2. General mathematical principles linking heteroplasmy variance to network dynamics.Wild-type and mutant copy numbers (A) and fused and unfused copy numbers (B) both move towards a line of steady states under a deterministic model, as indicated by arrows.In stochastic simulation, mean copy number (C) is initially slightly perturbed from the deterministic treatment of the system, and then remains constant, while mean heteroplasmy (D) remains invariant with time (see (Equation75)).In (E)-(H), we show that (Equation11) holds across many cellular circumstances: lines give analytic results, points are from stochastic simulation.Heteroplasmy variance behaviour is successfully predicted for varying mitophagy rate (E), steady state copy number (F), mutation sensing (G), and fusion rate (H).In (H), fusion and fission rates are redefined as γ → γ0M R and β → β0M where M and R denote the relative magnitude and ratio of the network rates, and γ0, β0 denote the nominal parameterizations of the fusion and fission rates respectively (see TableS1).FigureS2Fshows a sweep of M over the same logarithmic range when R = 1.See FigureS2H-P and TableS2for parameter sweeps numerically demonstrating the generality of the result for different mtDNA control modes.

Figure B1 .
Figure B1.Rate of de novo mutation accumulation is sensitive to the network state/mitophagy rate and copy number for a time-rescaled infinite sites Moran model.(A) An infinite sites Moran model where Q mutations occur per Moran step (see (Equation 12)).(B-D) Influence of our proposed intervention strategies.(B) Mean number of distinct mutations increases with the fraction of unfused mitochondria.This corresponds to a simple rescaling of time, so all but one of the parameterizations are shown in grey.(C) The mean number of mutations per mtDNA also increases with the fraction of unfused mitochondria.Inset shows that the mean number of mutations per mtDNA is independent of the number of mtDNAs per cell; values of n are the same as in (D).(D) Mean number of mutations per cell increases according to the population size of mtDNAs.Standard error in the mean is too small to visualise, so error bars are neglected, given 10 3 realizations.

Figure S1 .
Figure S1.Related to Method details.Phase portrait for an ODE representation of the system in (Equation15)-(Equation21). The system displays two steady states: a trivial steady state at s + f = 0, and a non-trivial steady state.At a steady state, both time derivatives in (Equation22) and (Equation23) vanish.Trajectories (blue) show the evolution of the system up to t = 1000.The direction and magnitude of the derivative at points in space are shown by red arrows.Trajectories can be seen either decaying to s = f = 0 or tending to infinity.γ = 0.01656, β = 33.12,ρ = 0.023, η = 1.1, α = 1.04.

Figure S2 .
Figure S2.Related to Figure 2, TableS2and Method details. (A) Deterministic dynamics of total copy number, which is controlled to a particular steady state value.(B) Defining a knock-down (KD) factor (k −1 = 0.1, 0.2 . . ., 1.0), the fission rate was rescaled to β → β/k (red) and the fusion rate to γ → γ/k (blue), causing a linear increase and decrease in total copy number respectively under a deterministic treatment.(C) Copy number variance for stochastic simulations initially increases, since all stochastic simulations begin with the same initial condition, but then plateaus since the steady state line is globally attracting.(D) Error in (Equation11) in a sweep over the feedback control strength, b.Dotted line denotes a 5% error according to (Equation99).(E) V(h) profile for the parameterization with the largest error in (D).(F) Sweeps of the network rate magnitude (see Figure2H).Heteroplasmy variance is approximately independent of absolute network rates over a broad range of network magnitudes.(G) Allowing fused species to be degraded with relative rate ξ (Equation82), stochastic simulations for heteroplasmy variance (markers) and (Equation83) (lines) are shown.Fused species degradation induces a linear correction to the heteroplasmy variance formula.(H-P) Error in (Equation11) (ansatz) and (Equation86) (ansatz indep network) for sweeps of the fusion and fission rates the corresponding feedback control functions in TableS2.Equations are accurate to at least 5% (blue regions) across large regions of parameter space, for many control laws.Fusion and fission rates are redefined as γ → γ0M R and β → β0M where M and R denote the magnitude and ratio of the network rates, and γ0, β0 denote the nominal parameterizations of the fusion and fission rates respectively.Summary statistics for 10 4 realizations, with initial condition h = 0.3 and evaluated at t = 500 days.Errors in V(h) (see (Equation99)) smaller than 5% are truncated and are shown as blue.Parameterizations where a deterministic steady state could not be found for an initial condition of h = 0.3 are shown in grey.Inset figures, where present, display the probability of fixation at h = 0.Where insets are not present, the probability of fixation is negligible.(H-M) When λ = λ(x) and P (h = 0) is low, (Equation11) performs well in the high-churn limit.(N-P) When µ = µ(x) and P (h = 0) is low, (Equation86) performs well in the high-churn limit.

Figure S4 .
Figure S4.Relating to Figure 3.A deterministic parameter sweep of fusion selectivity and the relative fusion rate for (Equation13), according to an ODE An ODE treatment allows smaller heteroplasmy changes to be probed without the need to resort to an infeasible number of stochastic simulations.Displaying the relative change in heteroplasmy (∆h) after t = 1000 days.We observe that a reduction in heteroplasmy is achieved at intermediate fusion rates at all non-zero fusion selectivities investigated.Grey denotes a change which is smaller than floating point precision.