Mitochondrial network state scales mtDNA genetic dynamics

Juvid Aryaman,∗,†,‡ Charlotte Bowles,§ Nick S. Jones,∗,∗∗,1 Iain G.
Johnston,††,§§,2

∗ Department of Mathematics, Imperial College London, London, SW7 2AZ, United
Kingdom

† Department of Clinical Neurosciences, University of Cambridge, Cambridge,
CB2 0QQ, United Kingdom

‡ MRC Mitochondrial Biology Unit, University of Cambridge, Cambridge, CB2 0XY,
United Kingdom

§ School of Biosciences, University of Birmingham, Birmingham, B15 2TT, United
Kingdom

∗∗ EPSRC Centre for the Mathematics of Precision Healthcare, Imperial College
London, London, SW7 2AZ, United Kingdom

†† Faculty of Mathematics and Natural Sciences, University of Bergen, 5007
Bergen, Norway

§§ Alan Turing Institute, London, NW1 2DB, United Kingdom

1<EMAIL_ADDRESS>
2<EMAIL_ADDRESS>

## Abstract

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 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 wide-reaching theoretical framework which
bridges mitochondrial network and genetic states. We show that, under a range
of conditions, the (genetic) rate of increase in heteroplasmy variance and de
novo mutation are proportionally modulated by the (physical) fraction of
unfused mitochondria, independently of the absolute fission-fusion 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 one of several leading hypotheses 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 (Wallace and Chalkia, 2013; Stewart
and Chinnery, 2015; Johnston, 2018; Aryaman et al., 2019). 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 (Aryaman et al., 2019).

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). Here, motivated
by a general picture of aging, we will largely focus on the setting of non-
dividing cells, which possess two mtDNA variants (although we will also
consider de novo mutation using simple statistical genetics models). The
variance of the distribution of heteroplasmies gives the fraction of cells
above a given pathological threshold (Figure 1B). Therefore heteroplasmy
variance is related to the number of dysfunctional cells above a phenotypic
threshold within a tissue, and both heteroplasmy mean and variance are
directly related to tissue physiology. Increases in heteroplasmy variance also
increase the number of cells below a given threshold heteroplasmy, which can
be advantageous in e.g. selecting low-heteroplasmy embryos in pre-implantation
genetic diagnosis for treating mitochondrial disease (Burgstaller et al.,
2014b; Johnston et al., 2015).

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). Heteroplasmy is depicted as an approximately normal distribution, as
this is the regime in which our approximations below hold: i.e. when the
probability of fixation is small. (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 ($\mathbb{V}(h)$)
approximately increases with time ($t$), mitophagy rate ($\mu$) and the
fraction of mitochondria that are unfused ($f_{s}$), and decreases with mtDNA
copy number ($n$). Importantly, $\mathbb{V}(h)$ does not depend on the
absolute magnitude of the fission-fusion rates. Also see Table S1 for a
summary of our key findings.

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) act in concert to
perform quality control on the mitochondrial population (Twig et al., 2008;
Johnston, 2018; Aryaman et al., 2019). Observations of pervasive intra-
mitochondrial mtDNA mutation (Morris et al., 2017) and universal heteroplasmy
in humans (Payne et al., 2012) suggest that the power of this quality control
may be limited. It has also been suggested that certain mtDNA mutations, such
as deletions (Kowald and Kirkwood, 2013, 2014, 2018) and some point mutations
(Samuels et al., 2013; Ye et al., 2014; Li et al., 2015; Lieber et al., 2019),
are under the influence of selective effects. However, genetic models without
selection have proven valuable in explaining the heteroplasmy dynamics both of
functional mutations (Elson et al., 2001; Taylor et al., 2003; Wonnapinij et
al., 2008) and polymorphisms without dramatic functional consequences (Birky
et al., 1983; Ye et al., 2014), and in common cases where mean heteroplasmy
shifts are small compared to changes in variances (for instance, in germline
development (Johnston et al., 2015) and post-mitotic tissues (Burgstaller et
al., 2014a)). Mean changes seem more likely in high-turnover tissues and when
mtDNA variants are genetically distant (Burgstaller et al., 2014a; Pan et al.,
2019), suggesting that neutral genetic theory may be useful in understanding
the dynamics of the set of functionally mild mutations which accumulate during
ageing. Neutral genetic theory also provides a valuable null model for
understanding mitochondrial genetic dynamics (Chinnery and Samuels, 1999;
Poovathingal et al., 2009; Johnston and Jones, 2016), potentially allowing us
to better understand and quantify when selection is present. There is thus a
set of open questions about how the physical dynamics of mitochondria affect
the genetic populations of mtDNA within and between cells under neutral
dynamics.

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)). These studies have suggested: that
clearance of damaged mtDNA can be assisted by high and funcitonally-selective
mitochondrial fusion, or by intermediate fusion and selective mitophagy (Mouli
et al., 2009); that physical transport of mitochondria can indirectly modulate
mitochondrial health through mitochondrial dynamics (Patel et al., 2013); that
fission-fusion dynamic rates modulate a trade-off between mutant proliferation
and removal (Tam et al., 2013, 2015); and that if fission is damaging,
decelerating fission-fusion cycles may improve mitochondrial quality (Figge et
al., 2012).

Despite providing valuable insights, these previous attempts to link
mitochondrial genetics and network dynamics, while important for breaking
ground, have centered around complex computer simulations, making it difficult
to deduce general laws and principles. Here, we address this lack of a general
theoretical framework linking mitochondrial dynamics and genetics. 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 (physical)
network state and (genetic) heteroplasmy statistics (Figure 1D). Our models
potentially allow rich interactions between mitochondrial genetic and network
dynamics, yet we find that a simple link emerges. For a broad range of
situations, 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. We discover that this result stems from the general
notion that fusion shields mtDNAs from turnover, since autophagy of large
fragments of the mitochondrial network are unlikely, and consequently rescales
time. Importantly, we used our model for network dynamics to 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)). Surprisingly, we find
the dependence of heteroplasmy statistics upon network state arises when the
mitochondrial population size is controlled through replication, and vanishes
when it is controlled through mitophagy, shedding new light on the
physiological importance of the mode of mtDNA control. We show 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
affect the rate of accumulation of de novo mutations (Johnston and Jones,
2016), see Table S1 for a summary of our key findings. We suggest that
pharmacological interventions which promote fusion, slow mitophagy and
increase copy number earlier in development may slow the rate of accumulation
of pathologically mutated cells, with implications for mitochondrial disease
and aging.

## Materials and Methods

### 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 seek to understand the influence of the mitochondrial
network upon mitochondrial genetics. The network state itself is influenced by
several factors including metabolic poise and the respiratory state of
mitochondria (Szabadkai et al., 2006; Hoitzing et al., 2015; Mishra and Chan,
2016), which we do not consider explicitly here. We consider the existence of
two mitochondrial alleles, wild-type ($W$) and mutant ($M$), existing within a
post-mitotic cell without cell division, with mtDNAs undergoing turnover (or
“relaxed replication” Stewart and Chinnery (2015)). 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

$\displaystyle X_{S}+X_{S}$ $\displaystyle\xrightarrow{\gamma}$ $\displaystyle
X_{F}+X_{F}$ (1) $\displaystyle X_{F}+X_{S}$
$\displaystyle\xrightarrow{\gamma}$ $\displaystyle X_{F}+X_{F}$ (2)
$\displaystyle X_{F}$ $\displaystyle\xrightarrow{\beta}$ $\displaystyle X_{S}$
(3)

where $X$ denotes either a wild-type ($W$) or a mutant ($M$) mtDNA (therefore
a set of chemical reactions analogous to Eq. (1)-(3) exist for both DNA
species). $\gamma$ and $\beta$ are the stochastic rate constants for fusion
and fission respectively.

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

$\displaystyle X_{S}$ $\displaystyle\xrightarrow{\lambda}$ $\displaystyle
X_{F}+X_{F}$ (4) $\displaystyle X_{F}$ $\displaystyle\xrightarrow{\lambda}$
$\displaystyle X_{F}+X_{F}$ (5) $\displaystyle X_{S}$
$\displaystyle\xrightarrow{\mu}$ $\displaystyle\emptyset$ (6)

where $\lambda$ and $\mu$ are the replication and mitophagy rates
respectively, which are shared by both $W$ and $M$ resulting in a so-called
‘neutral’ genetic model. Eq. (6) denotes removal of the species from the
system. The effect of allowing non-zero degradation of fused species is
discussed in Supporting Information (see Eq. (S68) and Figure S3E).
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, thereby
associating mitochondrial replication with fission (Lewis et al., 2016),
leaves our central result (Figure 1D) unchanged (see Supporting Information,
Eq. (S67)). 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:

$\displaystyle W_{F}+M_{S}$ $\displaystyle\xrightarrow{\gamma}$ $\displaystyle
W_{F}+M_{F}$ (7) $\displaystyle M_{F}+W_{S}$
$\displaystyle\xrightarrow{\gamma}$ $\displaystyle M_{F}+W_{F}$ (8)
$\displaystyle W_{S}+M_{S}$ $\displaystyle\xrightarrow{\gamma}$ $\displaystyle
W_{F}+M_{F}.$ (9)

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
Supporting Information, Rate renormalization for a discussion of rate
renormalization). We have neglected de novo mutation in the model description
above (although we will consider de novo mutation using a modified infinite
sites Moran model below).

We found that treating $\lambda=$ const led to instability in total copy
number (see Supporting Information, Constant rates yield unstable copy numbers
for a model describing mtDNA genetic and network dynamics), 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):

$\lambda=\lambda(w_{T},m_{T})=\mu+b(\kappa-(w_{T}+\delta m_{T}))$ (10)

where $w_{T}=w_{s}+w_{f}$ is the total wild-type copy number, and similarly
for $m_{T}$. The lower-case variables $w_{s}$, $w_{f}$, $m_{s}$, and $m_{f}$
denote the copy numbers of the corresponding chemical species ($W_{S}$,
$W_{F}$, $M_{S}$, and $M_{F}$). $b$ is a parameter which determines the
strength with which total copy number is controlled to a target copy number,
and $\kappa$ is a parameter which is indicative of (but not equivalent to) the
steady state copy number. $\delta$ 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\leq\delta<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. $\lambda$ may be re-
written as $\lambda=k_{1}+k_{2}w_{T}+k_{3}m_{T}$ for constants $k_{i}$, and so
only consists of 3 independent parameters. However we will retain $\lambda$ in
the form of Eq. (10) since the parameters $\mu$, $b$, $\kappa$, and $\delta$
have the distinct physiological meanings described above (Hoitzing et al.,
2017; Hoitzing, 2017). Furthermore, $\lambda$ may in general also depend on
other cellular features such as mitochondrial reactive oxygen species. Here,
we seek to explain mitochondrial behaviour under a simple set of governing
principles, but our approach can naturally be combined with a description of
these additional factors to build a more comprehensive model. 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 independent parameters
and 15 reactions, and captures the central property that mitochondria fragment
before degradation (Twig et al., 2008).

Throughout this work, we define heteroplasmy as the mutant allele fraction per
cell of a mitochondrially-encoded variant (Wonnapinij et al., 2008; Samuels et
al., 2010; Aryaman et al., 2019):

$h(\mathbf{x})=(m_{s}+m_{f})/(w_{s}+w_{f}+m_{s}+m_{f})$ (11)

where $\mathbf{x}=(w_{s},w_{f},m_{s},m_{f})$ is the state of the system (not
to be confused with mitochondrial “respiratory states”). Hence, a heteroplasmy
of $h=1$ denotes a cell with 100% mutant mtDNA (i.e. a homoplasmic cell in the
mutant allele). Arguably, “mutant allele fraction” would be a more precise
description of Eq.(11) but we retain the use of heteroplasmy for consistency.
To convert to a definition of heteroplasmy which is maximal when the mutant
allele fraction is 50%, one may simply use the conversion
$0.5-|h(\mathbf{x})-0.5|$.

#### Statistical Analysis

In Figures S3B, S4A-I, we compare Eq. (13) and Eq. (S72) to stochastic
simulations, for various parametrizations and replication/degradation rates.
To quantify the accuracy of these equations in predicting $\mathbb{V}(h,t)$,
we define the following error metric $\epsilon$

$\epsilon=\left|1-\frac{\dot{\mathbb{V}}(h,t)_{\text{Th}}}{\mathbb{E}_{t}(\dot{\mathbb{V}}(h,t)_{\text{Sim}})}\right|$
(12)

where $\dot{\mathbb{V}}(h,t)$ is the time derivative of heteroplasmy variance
with subscripts denoting theory (Th) and simulation (Sim). An expectation over
time ($\mathbb{E}_{t}$) is taken for the stochastic simulations, whereas
$\dot{\mathbb{V}}(h,t)$ is a scalar quantity for Eq. (13) and Eq. (S72).

### Data Availability

Code for simulations and analysis can be accessed at
https://GitHub.com/ImperialCollegeLondon/MitoNetworksGenetics

## Results

### 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 Eqs.
(1)–(10), by converting the reactions into an analogous set of four coupled
ordinary differential equations (see Eqs. (S29)–(S32)), and choosing a
biologically-motivated approximate parametrization (which we will term the
‘nominal’ parametrization, see Supporting Information, Choice of nominal
parametrization, and Table S2). Figures 2A-B show that copy numbers of each
individual species change in time such that the state approaches a line of
steady states (Eqs. (S34)–(S36)), 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 Eq. (10) which decreases with total copy number when $w_{T}+\delta
m_{T}>\kappa$ and vice versa, guiding the total population to a fixed total
copy number. Varying the fission ($\beta$) and fusion ($\gamma$) 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 Eqs. (1)–(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
$\delta=1$ for our nominal parametrization (i.e. the line of steady states is
also a line of constant copy number). Choosing $\delta\neq 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 S3A). 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.

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 Eq. (S61)). In (E)-(H), we show that Eq.
(13) 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 $\gamma\rightarrow\gamma_{0}MR$ and
$\beta\rightarrow\beta_{0}M$ where $M$ and $R$ denote the relative magnitude
and ratio of the network rates, and $\gamma_{0},\beta_{0}$ denote the nominal
parametrizations of the fusion and fission rates respectively (see Table S2).
Figure S3D shows a sweep of $M$ over the same logarithmic range when $R=1$.
See Figure S4A-I and Table S3 for parameter sweeps numerically demonstrating
the generality of the result for different mtDNA control modes.

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.

To this end, we analytically explored the influence of mitochondrial dynamics
on mtDNA variability. Assuming that the state of the system above is
initialised at its deterministic steady state
($\mathbf{x}(t=0)=\mathbf{x}_{\text{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 Supporting
Information). 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 Eq.
(S63), which may be approximated for sufficiently short times as

$\mathbb{V}(h)\approx\left.f_{s}(\mathbf{x})\frac{2\mu
t}{n(\mathbf{x})}h(\mathbf{x})(1-h(\mathbf{x}))\right|_{\mathbf{x}=\mathbf{x}_{\text{ss}}}.$
(13)

Here, $\mathbb{V}(h)$ is the variance of heteroplasmy, $\mu$ is the mitophagy
rate, $n(\mathbf{x})$ is the total copy number and $f_{s}(\mathbf{x})$ is the
fraction of unfused (singleton) mtDNAs, and is thus a measure of the
fragmentation of the mitochondrial network. $\mathbf{x}_{\text{ss}}$ is the
(deterministic) steady state of the system. Eq. (13) demonstrates that mtDNA
heteroplasmy variance increases approximately linearly with time ($t$) at a
rate scaled by the fraction of unfused mitochondria, mitophagy rate, and
inverse population size. We find that Eq. (13) 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, Figure S5).

To our knowledge, Eq. (13) 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 $\mathbb{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 $\mu$, 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 $\mathbb{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 below) 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
Eq. (13) is able to recapitulate the rate of heteroplasmy variance increase
across different values of $\delta$, which are hypothesized to correspond to
different replicative sensing strengths of different mitochondrial mutations
(Hoitzing, 2017). We also show in Figures S3B&C that Eq. (13) is robust to the
choice of feedback control strength $b$ in Eq. (10). $n(\mathbf{x})$,
$f(\mathbf{x})$, and $h(\mathbf{x})$ in Eq.(13) are not independent degrees of
freedom in this model: they are functions of the state vector $\mathbf{x}$,
where $\mathbf{x}$ is determined by the parametrization and initial conditions
of the model. Hence, the parameter sweeps in Figure 2E-H and Figures S3B&C
also implicitly vary over these functions of state by varying the steady state
$\mathbf{x}_{\text{ss}}$.

In Eq. (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 S3D). This observation,
which contrasts with the model of (Tam et al., 2013, 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
$\beta/\mu\approx 10^{3}$ (see Table S2), resulting in the existence of a
separation of timescales between network and genetic processes. In the
derivation of Eq. (13), we have assumed that fission-fusion rates are
infinite, which simplifies $\mathbb{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, Eq. (13) breaks down and $\mathbb{V}(h)$ gains
dependence upon the magnitude of these rates (see Figure S4A). This regime is,
however, for network rates which are approximately 100 times smaller than the
biologically-motivated nominal parametrization 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.

Eq. (13) can be viewed as describing the “quasi-stationary state” where the
probability of extinction of either allele is negligible (Johnston and Jones,
2016). On longer timescales, or if mtDNA half-life is short (Poovathingal et
al., 2012), the probability of fixation becomes appreciable. In this case, Eq.
(13) over-estimates $\mathbf{V}(h)$ as heteroplasmy variance gradually becomes
sub-linear with time, see Figure S5C&D. This is evident through inspection of
Eq. (S63), which shows that cellular trajectories which reach $h=0$ or $h=1$
cease to diffuse in heteroplasmy space, and so heteroplasmy variance cannot
increase indefinitely. Consequently, the depiction of heteroplasmy variance in
Fig. 1B,D as being approximately normally distributed corresponds to the
regime in which our approximation holds, and is a valid subset of the
behaviours displayed by heteroplasmy dynamics under more sophisticated models
(e.g. the Kimura distribution Kimura (1955); Wonnapinij et al. (2008)).
Further analytical developments may be possible to take into account
extinction (e.g. see Wonnapinij et al. (2008); Assaf and Meerson (2010)).
However, the linear regime for heteroplasmy variance has been observed to be a
substantial component of mtDNA dynamics in e.g. mouse oocytes (Burgstaller et
al., 2018).

#### 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.
$\lambda=\lambda(\mathbf{x})$, $\mu=$ const), when the fusion and fission
rates were high and the fixation probability ($P(h=0)$ or $P(h=1)$) was
negligible, Eq. (13) accurately described $\mathbb{V}(h)$ across all of the
replication rates investigated, see Figure S4A-F. The same mathematical
argument to show Eq. (13) for the replication rate in Eq. (10) may be applied
to these alternative replication rates where a closed-form solution for the
deterministic steady state may be written down (see Supporting Information,
Deriving an ODE description of the mitochondrial network system).
Interestingly, when copy number is controlled through the degradation rate
(i.e. $\lambda=$const, $\mu=\mu(\mathbf{x})$), heteroplasmy variance loses its
dependence upon network state entirely and the $f_{s}$ term is lost from Eq.
(13) (see Eq. (S72) and Figure S4G-I). A similar mathematical argument was
applied to reveal how this dependence is lost (see Supporting Information,
Proof of heteroplasmy relation for linear feedback control).

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 Supporting Information, A modified Moran
process may account for the alternative forms of heteroplasmy variance
dynamics under different models of genetic mtDNA control). 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 Eq. (S84) and surrounding discussion).

### 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 ($\mu$), decreases with
total mtDNA copy number per cell ($n$) and increases with the fraction of
unfused mitochondria (termed “singletons”, $f_{s}$), see Eq. (13). 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 against mutant expansions

In order to explore the role of the mitochondrial network in the accumulation
of de novo mutations, we invoked an infinite sites Moran model (Kimura, 1969)
(see Figure 3A). 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

$Q\sim\operatorname{\text{Binomial}}(L_{\text{mtDNA}},\eta)$ (14)

where $\operatorname{\text{Binomial}}(N,p)$ is a binomial random variable with
$N$ trials and probability $p$ of success. $L_{\text{mtDNA}}=16569$ is the
length of mtDNA in base pairs and $\eta=5.6\times 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. In
Supporting Information, Eq. (S83), we argue that when population size is
controlled in the replication rate, the inter-event rate ($\Gamma$) of the
Moran process is effectively rescaled by the fraction of unfused mitochondria,
i.e. $\Gamma=\mu nf_{s}$, which we apply here.

Figure 3: 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 Eq. (14)). (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 parametrizations 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 3B shows that in the infinite sites model, the consequence of Eq. (S83)
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 3C). 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 3B&C therefore correspond to
a rescaling of time i.e. stretching of the time-axis. The absolute number of
mutations predicted in Figure 3B may over-estimate the true number of
mutations per cell (and of course depends on our choice of mutation rate),
since a subset of mutations will experience either positive or negative
selection. However, quantification of the number of distinct mitochondrial
mutants in single cells remains under-explored, as most mutations will have a
variant allele fraction close to 0% or 100% (Birky et al., 1983), which are
challenging to measure, especially through bulk sequencing.

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. Although knock-out studies present difficulties in extending their
insights into the physiological case, the authors observed that fragmentation
of the mitochondrial network induced severe depletion of mtDNA copy number
(which we also observed in Figure S2B). Furthermore, the authors observed that
the number of mutations per base pair increased upon fragmentation, which we
also observed in the infinite sites model where fragmentation effectively
results in a faster turnover of mtDNA (Figure 3C).

Our models predict that promoting mitochondrial fusion has a two-fold effect:
firstly, it slows the increase of heteroplasmy variance (see Eq. (13) and
Figure 2H); secondly, it reduces the rate of accumulation of distinct
mutations (see Figure 3B&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, this simple
model suggests that promoting fusion earlier in development (assuming mean
heteroplasmy is low) could slow down the accumulation and spread of
mitochondrial mutations, and perhaps slow aging.

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 (see below), 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 (see below).

#### Targeting mitophagy rate against mutant expansions

Alterations in the mitophagy rate $\mu$ have a comparable effect to changes in
$f_{s}$ in terms of reducing the rate of heteroplasmy variance (see Eq. (13))
and the rate of de novo mutation (Figure 3B&C) since they both serve to
rescale time. Our theory therefore suggests that inhibition of basal mitophagy
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 against mutant expansions

To investigate the role of mtDNA copy number (mtCN) on the accumulation of de
novo mutations, we set $f_{s}=1$ such that $\Gamma=\mu 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 3C, inset). However, as the population
size becomes larger, the total number of distinct mutations increases
accordingly (Figure 3D). 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 Eq.
(13)). 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 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 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 ($\gamma$) in the mutant
analogues of Eq. (1),(2) and Eqs. (7)–(9) by writing

$\gamma\rightarrow\gamma/(1+\epsilon_{f})$ (15)

for all fusion reactions involving 1 or more mutant mitochondria where
$\epsilon_{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

$M_{s}\xrightarrow{\mu(1+\epsilon_{m})}\emptyset$ (16)

for $\epsilon_{m}>0$.

In these two settings, we explore how varying the fusion rate for a given
selectivity ($\epsilon_{f}$ and $\epsilon_{m}$) affects the extent of
reduction in mean heteroplasmy. Figure 4A shows that, in the context of
selective fusion ($\epsilon_{f}>0$) and non-selective mitophagy
($\epsilon_{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 S7) 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.

Figure 4: 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 $\epsilon_{f}=\epsilon_{m}=0$, for $10^{4}$
iterations. (A) For selective fusion (see Eq. (15)), for each value of fusion
selectivity ($\epsilon_{f}$), the fusion rate ($\gamma$) was varied relative
to the nominal parametrization (see Table S2). When $\epsilon_{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 S7). (B) For selective mitophagy (see
Eq. (16)), when mitophagy selectivity $\epsilon_{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.

In contrast, in Figure 4B, in the context of non-selective fusion
($\epsilon_{f}=0$) and selective mitophagy ($\epsilon_{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 physiology – the 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. 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.

Our bottom-up modeling approach allows for potentially complex interactions
between the physical (network) and genetic mitochondrial states of the cell,
yet a simple connection emerged from our analysis. We found, for a wide class
of models of post-mitotic cells, that the rate of linear increase of
heteroplasmy variance is modulated in proportion to the fraction of unfused
mitochondria (see Eq. (13)). The general notion that mitochondrial fusion
shields mtDNAs from turnover, and consequently serves to rescale time, emerges
from our analysis. This rescaling of time only holds when mitochondrial copy
numbers are controlled through a state-dependent replication rate, and
vanishes if copy numbers are controlled through a state-dependent mitophagy
rate. 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 3B&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 observations provide 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, through building a model that incorporates
mitochondrial dynamics, 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
information about heteroplasmy dynamics via the mitochondrial network, without
the need to quantify absolute fission-fusion rates (for instance through
confocal micrographs to quantify the fraction of unfused mitochondria). 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 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).

We have, broadly speaking, considered 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 Eq. (15) and Eq. (16)).
If mitophagy is able to change mean heteroplasmy, 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 paid particular attention to the case of post-mitotic tissues, since
these tissues are important for understanding the role of mitochondrial
mutations in healthy aging (Khrapko and Vijg, 2009; Kauppila et al., 2017). A
typical rate of increase of heteroplasmy variance predicted by Eq.(13) given
our nominal parametrization (Table S2) is
$V^{\prime}(h)/t=\mathbb{V}(h)/(\mathbb{E}(h)(1-\mathbb{E}(h))t)=2\mu
f_{s}/n\approx 2.3\times 10^{-5}$ day-1 ($f_{s}=0.5$, $n=1000$). This value
accounts for the accumulation of heteroplasmy variance which is attributable
to turnover of the mitochondrial population in a post-mitotic cell. However,
in the most general case, cell division is also able to induce substantial
heteroplasmy variance. For example, $V^{\prime}(h)/t$ has been measured in
model organism germlines to be approximately $9\times 10^{-4}$ day-1 in
Drosophila (Solignac et al., 1987; Johnston and Jones, 2016), $9\times
10^{-4}$ day-1 in NZB/BALB mice (Wai et al., 2008; Wonnapinij et al., 2008;
Johnston and Jones, 2016), and $2\times 10^{-4}$ day-1 in single LE and HB
mouse oocytes (Burgstaller et al., 2018). We see that these rates of increase
in heteroplasmy variance are approximately an order of magnitude larger than
predictions from our model of purely quiescent turnover, given our nominal
parametrisation. Whilst larger mitophagy rates may also potentially induce
larger values for $V^{\prime}(h)/t$ (see Poovathingal et al. (2012), and
Figure S5C, corrsponding to $V^{\prime}(h)/t\approx 3.5\times 10^{-4}$ day-1)
it is clear that partitioning noise (or “vegetative segregation”, Stewart and
Chinnery (2015)) is also an important source of variance in heteroplasmy
dynamics (Johnston et al., 2015). Quantification of heteroplasmy variance in
quiescent tissues remains an under-explored area, despite its importance in
understanding healthy ageing (Kauppila et al., 2017; Aryaman et al., 2019).

Our findings reveal some apparent differences with previous studies which link
mitochondrial genetics with network dynamics (see Table S4). Firstly, Tam et
al. (2013, 2015) found that slower fission-fusion dynamics resulted in larger
increases in heteroplasmy variance with time, in contrast to Eq. (13) which
only depends on fragmentation state and not absolute network rates. The
simulation approach of Tam et al. (2013, 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., 2014a). 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, 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. (2013, 2015) is the inclusion of cellular
subcompartments which induces spatial effects which we do not consider here.
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 4A). 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 highest 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 4A.
Furthermore, Mouli et al. (2009) also found that when fusion is non-selective
and mitophagy is selective, intermediate fusion rates are optimal whereas
Figure 4B 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 4B 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. 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. (2013, 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}$, $\mu$ and $n$, with
heteroplasmy quantification, would allow verification of Eq. (13) 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}$, $\mu$ and $n$, followed by e.g. destructive
single-cell whole-genome sequencing of mtDNA would allow validation of how
$\mu$, $n$ and $f_{s}$ influence $\mathbb{V}(h)$ and the rate of de novo
mutation (see Figure 3). 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.

### Acknowledgments

We would like to thank Hanne Hoitzing, Thomas McGrath, Abhishek Deshpande, and
Ferdinando Insalata for useful discussions. Simulations were performed using
the Imperial College High Performance Computing Service. J.A. acknowledges
grant support from the BBSRC (BB/J014575/1), and the MRC Mitochondrial Biology
Unit (MC_UP_1501/2). N.S.J. acknowledges grant support from the BHF
(RE/13/2/30182) and EPSRC (EP/N014529/1). I.G.J. acknowledges funding from ERC
StG 805046 (EvoConBiO) and a Turing Fellowship from the Alan Turing Institute.

## References

  * Acton et al. (2007) Acton, B., I. Lai, X. Shang, A. Jurisicova, and R. Casper, 2007 Neutral mitochondrial heteroplasmy alters physiological function in mice. Biol. Reprod. 77: 569–576.
  * Aryaman et al. (2017) Aryaman, J., I. G. Johnston, and N. S. Jones, 2017 Mitochondrial DNA density homeostasis accounts for a threshold effect in a cybrid model of a human mitochondrial disease. Biochem. J. 474: 4019–4034.
  * Aryaman et al. (2019) Aryaman, J., I. G. Johnston, and N. S. Jones, 2019 Mitochondrial heterogeneity. Front. Genet. 9: 718.
  * Assaf and Meerson (2010) Assaf, M. and B. Meerson, 2010 Extinction of metastable stochastic populations. Phys. Rev. E 81: 021116.
  * Bagwan et al. (2018) Bagwan, N., E. Bonzon-Kulichenko, E. Calvo, A. V. Lechuga-Vieco, S. Michalakopoulos, et al., 2018 Comprehensive quantification of the modified proteome reveals oxidative heart damage in mitochondrial heteroplasmy. Cell Rep. 23: 3685–3697.
  * Birky et al. (1983) Birky, C. W., T. Maruyama, and P. Fuerst, 1983 An approach to population and evolutionary genetic theory for genes in mitochondria and chloroplasts, and some results. Genetics 103: 513–527.
  * Boudoures et al. (2017) Boudoures, A. L., J. Saben, A. Drury, S. Scheaffer, Z. Modi, et al., 2017 Obesity-exposed oocytes accumulate and transmit damaged mitochondria due to an inability to activate mitophagy. Dev. Biol. 426: 126–138.
  * Burgstaller et al. (2014a) Burgstaller, J. P., I. G. Johnston, N. S. Jones, J. Albrechtová, T. Kolbe, et al., 2014a MtDNA segregation in heteroplasmic tissues is common in vivo and modulated by haplotype differences and developmental stage. Cell Rep. 7: 2031–2041.
  * Burgstaller et al. (2014b) Burgstaller, J. P., I. G. Johnston, and J. Poulton, 2014b Mitochondrial DNA disease and developmental implications for reproductive strategies. Mol. Hum. Reprod. 21: 11–22.
  * Burgstaller et al. (2018) Burgstaller, J. P., T. Kolbe, V. Havlicek, S. Hembach, J. Poulton, et al., 2018 Large-scale genetic analysis reveals mammalian mtDNA heteroplasmy dynamics and variance increase through lifetimes and generations. Nat. Commun. 9: 2488.
  * Cagalinec et al. (2013) Cagalinec, M., D. Safiulina, M. Liiv, J. Liiv, V. Choubey, et al., 2013 Principles of the mitochondrial fusion and fission cycle in neurons. J. Cell Sci. 126: 2187–2197.
  * Capps et al. (2003) Capps, G. J., D. C. Samuels, and P. F. Chinnery, 2003 A model of the nuclear control of mitochondrial DNA replication. J. Theor. Biol. 221: 565–583.
  * Chen et al. (2010) Chen, H., M. Vermulst, Y. E. Wang, A. Chomyn, T. A. Prolla, et al., 2010 Mitochondrial fusion is required for mtDNA stability in skeletal muscle and tolerance of mtDNA mutations. Cell 141: 280–289.
  * Chinnery and Samuels (1999) Chinnery, P. F. and D. C. Samuels, 1999 Relaxed replication of mtDNA: a model with implications for the expression of disease. Am. J. Hum. Genet. 64: 1158.
  * Constable et al. (2016) Constable, G. W., T. Rogers, A. J. McKane, and C. E. Tarnita, 2016 Demographic noise can reverse the direction of deterministic selection. Proc. Natl. Acad. Sci. USA p. 201603693.
  * Cree et al. (2008) Cree, L., S. Patel, A. Pyle, S. Lynn, D. Turnbull, et al., 2008 Age-related decline in mitochondrial DNA copy number in isolated human pancreatic islets. Diabetologia 51: 1440–1443.
  * Dai et al. (2013) Dai, Y., K. Zheng, J. Clark, R. H. Swerdlow, S. M. Pulst, et al., 2013 Rapamycin drives selection against a pathogenic heteroplasmic mitochondrial DNA mutation. Hum. Mol. Genet. 23: 637–647.
  * de Almeida et al. (2017) de Almeida, M. J., L. L. Luchsinger, D. J. Corrigan, L. J. Williams, and H.-W. Snoeck, 2017 Dye-independent methods reveal elevated mitochondrial mass in hematopoietic stem cells. Cell Stem Cell 21: 725–729.
  * Dölle et al. (2016) Dölle, C., I. Flønes, G. S. Nido, H. Miletic, N. Osuagwu, et al., 2016 Defective mitochondrial DNA homeostasis in the substantia nigra in Parkinson disease. Nat. Commun. 7: 13548.
  * Durham et al. (2007) Durham, S. E., D. C. Samuels, L. M. Cree, and P. F. Chinnery, 2007 Normal levels of wild-type mitochondrial DNA maintain cytochrome c oxidase activity for two pathogenic mitochondrial DNA mutations but not for m.3243A$\rightarrow$G. Am. J. Hum. Genet. 81: 189–195.
  * Elson et al. (2001) Elson, J., D. Samuels, D. Turnbull, and P. Chinnery, 2001 Random intracellular drift explains the clonal expansion of mitochondrial dna mutations with age. Am. J. Hum. Genet. 68: 802–806.
  * Figge et al. (2012) Figge, M. T., A. S. Reichert, M. Meyer-Hermann, and H. D. Osiewacz, 2012 Deceleration of fusion–fission cycles improves mitochondrial quality control during aging. PLoS Comp. Biol. 8: e1002576.
  * Gardiner (1985) Gardiner, C., 1985 Stochastic Methods: A Handbook for the Natural and Social Sciences. Springer Series in Synergetics, Springer Berlin Heidelberg.
  * Gillespie (1976) Gillespie, D. T., 1976 A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comp. Phys. 22: 403–434.
  * Gillespie (2007) Gillespie, D. T., 2007 Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem. 58: 35–55.
  * Grima (2010) Grima, R., 2010 An effective rate equation approach to reaction kinetics in small volumes: Theory and application to biochemical reactions in nonequilibrium steady-state conditions. J. Chem. Phys. 133: 07B604.
  * Grima et al. (2011) Grima, R., P. Thomas, and A. V. Straube, 2011 How accurate are the nonlinear chemical Fokker-Planck and chemical Langevin equations? J. Chem. Phys. 135: 084103.
  * Hoitzing (2017) Hoitzing, H., 2017 Controlling mitochondrial dynamics: population genetics and networks. Ph.D. thesis, Imperial College London.
  * Hoitzing et al. (2017) Hoitzing, H., P. A. Gammage, M. Minczuk, I. G. Johnston, and N. S. Jones, 2017 Energetic costs of cellular and therapeutic control of stochastic mtDNA populations. bioRxiv p. 145292.
  * Hoitzing et al. (2015) Hoitzing, H., I. G. Johnston, and N. S. Jones, 2015 What is the function of mitochondrial networks? A theoretical assessment of hypotheses and proposal for future research. BioEssays 37: 687–700.
  * Jacobs (2010) Jacobs, K., 2010 Stochastic processes for physicists: understanding noisy systems. Cambridge University Press.
  * Johnston (2018) Johnston, I. G., 2018 Tension and resolution: Dynamic, evolving populations of organelle genomes within plant cells. Mol. Plant .
  * Johnston et al. (2015) Johnston, I. G., J. P. Burgstaller, V. Havlicek, T. Kolbe, T. Rülicke, et al., 2015 Stochastic modelling, bayesian inference, and new in vivo measurements elucidate the debated mtDNA bottleneck mechanism. eLife 4: e07464.
  * Johnston and Jones (2016) Johnston, I. G. and N. S. Jones, 2016 Evolution of cell-to-cell variability in stochastic, controlled, heteroplasmic mtDNA populations. Am. J. Hum. Genet. 99: 1150–1162.
  * Kandul et al. (2016) Kandul, N. P., T. Zhang, B. A. Hay, and M. Guo, 2016 Selective removal of deletion-bearing mitochondrial DNA in heteroplasmic drosophila. Nat. Commun. 7: 13100.
  * Kauppila et al. (2017) Kauppila, T. E., J. H. Kauppila, and N.-G. Larsson, 2017 Mammalian mitochondria and aging: an update. Cell Metab. 25: 57–71.
  * Khrapko and Vijg (2009) Khrapko, K. and J. Vijg, 2009 Mitochondrial DNA mutations and aging: devils in the details? Trends Genet. 25: 91–98.
  * Kim et al. (2007) Kim, I., S. Rodriguez-Enriquez, and J. J. Lemasters, 2007 Selective degradation of mitochondria by mitophagy. Arch. Biochem. Biophys. 462: 245–253.
  * Kimura (1955) Kimura, M., 1955 Solution of a process of random genetic drift with a continuous model. Proc. Natl. Acad. Sci. USA 41: 144.
  * Kimura (1969) Kimura, M., 1969 The number of heterozygous nucleotide sites maintained in a finite population due to steady flux of mutations. Genetics 61: 893–903.
  * Kowald and Kirkwood (2018) Kowald, A. and T. Kirkwood, 2018 Resolving the enigma of the clonal expansion of mtdna deletions. Genes 9: 126.
  * Kowald and Kirkwood (2013) Kowald, A. and T. B. Kirkwood, 2013 Mitochondrial mutations and aging: random drift is insufficient to explain the accumulation of mitochondrial deletion mutants in short-lived animals. Aging Cell 12: 728–731.
  * Kowald and Kirkwood (2014) Kowald, A. and T. B. Kirkwood, 2014 Transcription could be the key to the selection advantage of mitochondrial deletion mutants in aging. Proc. Natl. Acad. Sci. USA 111: 2972–2977.
  * Kowald and Klipp (2014) Kowald, A. and E. Klipp, 2014 Mathematical models of mitochondrial aging and dynamics. In Prog. Mol. Biol. Transl. Sci., volume 127, pp. 63–92, Elsevier.
  * Kukat et al. (2011) Kukat, C., C. A. Wurm, H. Spåhr, M. Falkenberg, N.-G. Larsson, et al., 2011 Super-resolution microscopy reveals that mammalian mitochondrial nucleoids have a uniform size and frequently contain a single copy of mtDNA. Proc. Natl. Acad. Sci. USA 108: 13534–13539.
  * Lewis et al. (2016) Lewis, S. C., L. F. Uchiyama, and J. Nunnari, 2016 ER-mitochondria contacts couple mtDNA synthesis with mitochondrial division in human cells. Science 353: aaf5549.
  * Li et al. (2015) Li, M., R. Schröder, S. Ni, B. Madea, and M. Stoneking, 2015 Extensive tissue-related and allele-related mtDNA heteroplasmy suggests positive selection for somatic mutations. Proc. Natl. Acad. Sci. USA 112: 2491–2496.
  * Lieber et al. (2019) Lieber, T., S. P. Jeedigunta, J. M. Palozzi, R. Lehmann, and T. R. Hurd, 2019 Mitochondrial fragmentation drives selective removal of deleterious mtDNA in the germline. Nature p. 1.
  * López-Otín et al. (2013) López-Otín, C., M. A. Blasco, L. Partridge, M. Serrano, and G. Kroemer, 2013 The hallmarks of aging. Cell 153: 1194–1217.
  * McWilliams et al. (2016) McWilliams, T. G., A. R. Prescott, G. F. Allen, J. Tamjar, M. J. Munson, et al., 2016 mito-QC illuminates mitophagy and mitochondrial architecture in vivo. J. Cell Biol. pp. jcb–201603039.
  * McWilliams et al. (2018) McWilliams, T. G., A. R. Prescott, L. Montava-Garriga, G. Ball, F. Singh, et al., 2018 Basal mitophagy occurs independently of PINK1 in mouse tissues of high metabolic demand. Cell Metab. 27: 439–449.
  * Medeiros et al. (2018) Medeiros, T. C., R. L. Thomas, R. Ghillebert, and M. Graef, 2018 Autophagy balances mtDNA synthesis and degradation by DNA polymerase POLG during starvation. J. Cell Biol. pp. jcb–201801168.
  * Mengel-From et al. (2014) Mengel-From, J., M. Thinggaard, C. Dalgård, K. O. Kyvik, K. Christensen, et al., 2014 Mitochondrial DNA copy number in peripheral blood cells declines with age and is associated with general health among elderly. Hum. Genet. 133: 1149–1159.
  * Mishra and Chan (2016) Mishra, P. and D. C. Chan, 2016 Metabolic regulation of mitochondrial dynamics. J. Cell Biol. 212: 379–387.
  * Morris et al. (2017) Morris, J., Y.-J. Na, H. Zhu, J.-H. Lee, H. Giang, et al., 2017 Pervasive within-mitochondrion single-nucleotide variant heteroplasmy as revealed by single-mitochondrion sequencing. Cell Rep. 21: 2706–2713.
  * Mouli et al. (2009) Mouli, P. K., G. Twig, and O. S. Shirihai, 2009 Frequency and selectivity of mitochondrial fusion are key to its quality maintenance function. Biophys. J. 96: 3509–3518.
  * Narendra et al. (2008) Narendra, D., A. Tanaka, D.-F. Suen, and R. J. Youle, 2008 Parkin is recruited selectively to impaired mitochondria and promotes their autophagy. J. Cell Biol. 183: 795–803.
  * Pan et al. (2019) Pan, J., L. Wang, C. Lu, Y. Zhu, Z. Min, et al., 2019 Matching mitochondrial DNA haplotypes for circumventing tissue-specific segregation bias. iScience 13: 371–379.
  * Parsons and Rogers (2017) Parsons, T. L. and T. Rogers, 2017 Dimension reduction for stochastic dynamical systems forced onto a manifold by large drift: a constructive approach with examples from theoretical biology. J. Phys. A 50: 415601.
  * Patel et al. (2013) Patel, P. K., O. Shirihai, and K. C. Huang, 2013 Optimal dynamics for quality control in spatially distributed mitochondrial networks. PLoS Comp. Biol. 9: e1003108.
  * Payne et al. (2012) Payne, B. A., I. J. Wilson, P. Yu-Wai-Man, J. Coxhead, D. Deehan, et al., 2012 Universal heteroplasmy of human mitochondrial DNA. Hum. Mol. Genet. 22: 384–390.
  * Perier et al. (2013) Perier, C., A. Bender, E. García-Arumí, M. J. Melia, J. Bove, et al., 2013 Accumulation of mitochondrial DNA deletions within dopaminergic neurons triggers neuroprotective mechanisms. Brain 136: 2369–2378.
  * Picard et al. (2014) Picard, M., J. Zhang, S. Hancock, O. Derbeneva, R. Golhar, et al., 2014 Progressive increase in mtDNA 3243A$>$G heteroplasmy causes abrupt transcriptional reprogramming. Proc. Natl. Acad. Sci. USA 111: E4033–E4042.
  * Poovathingal et al. (2009) Poovathingal, S. K., J. Gruber, B. Halliwell, and R. Gunawan, 2009 Stochastic drift in mitochondrial DNA point mutations: a novel perspective ex silico. PLoS Comput. Biol. 5: e1000572.
  * Poovathingal et al. (2012) Poovathingal, S. K., J. Gruber, L. Lakshmanan, B. Halliwell, and R. Gunawan, 2012 Is mitochondrial DNA turnover slower than commonly assumed? Biogerontology 13: 557–564.
  * Rossignol et al. (2003) Rossignol, R., B. Faustin, C. Rocher, M. Malgat, J. Mazat, et al., 2003 Mitochondrial threshold effects. Biochem. J. 370: 751–762.
  * Salk et al. (2018) Salk, J. J., M. W. Schmitt, and L. A. Loeb, 2018 Enhancing the accuracy of next-generation sequencing for detecting rare and subclonal mutations. Nature Rev. Genet. .
  * Samuels et al. (2013) Samuels, D. C., C. Li, B. Li, Z. Song, E. Torstenson, et al., 2013 Recurrent tissue-specific mtDNA mutations are common in humans. PLoS Genet. 9: e1003929.
  * Samuels et al. (2010) Samuels, D. C., P. Wonnapinij, L. Cree, and P. F. Chinnery, 2010 Reassessing evidence for a postnatal mitochondrial genetic bottleneck. Nat. Genet. 42: 471–472.
  * Schon et al. (2012) Schon, E. A., S. DiMauro, and M. Hirano, 2012 Human mitochondrial DNA: roles of inherited and somatic mutations. Nat. Rev. Genet. 13: 878–890.
  * Sharpley et al. (2012) Sharpley, M. S., C. Marciniak, K. Eckel-Mahan, M. McManus, M. Crimi, et al., 2012 Heteroplasmy of mouse mtDNA is genetically unstable and results in altered behavior and cognition. Cell 151: 333–343.
  * Solignac et al. (1987) Solignac, M., J. Génermont, M. Monnerot, and J.-C. Mounolou, 1987 Drosophila mitochondrial genetics: evolution of heteroplasmy through germ line cell divisions. Genetics 117: 687–696.
  * Stewart and Chinnery (2015) Stewart, J. B. and P. F. Chinnery, 2015 The dynamics of mitochondrial DNA heteroplasmy: implications for human health and disease. Nature Rev. Genet. 16: 530–542.
  * Stewart and Larsson (2014) Stewart, J. B. and N.-G. Larsson, 2014 Keeping mtDNA in shape between generations. PLoS Genet. 10: e1004670.
  * Suen et al. (2010) Suen, D.-F., D. P. Narendra, A. Tanaka, G. Manfredi, and R. J. Youle, 2010 Parkin overexpression selects against a deleterious mtDNA mutation in heteroplasmic cybrid cells. Proc. Natl. Acad. Sci. USA 107: 11835–11840.
  * Sukhorukov et al. (2012) Sukhorukov, V. M., D. Dikov, A. S. Reichert, and M. Meyer-Hermann, 2012 Emergence of the mitochondrial reticulum from fission and fusion dynamics. PLoS Comput. Biol. 8: e1002745.
  * Szabadkai et al. (2006) Szabadkai, G., A. M. Simoni, K. Bianchi, D. De Stefani, S. Leo, et al., 2006 Mitochondrial dynamics and Ca2+ signaling. Biochim. Biophys. Acta. 1763: 442–449.
  * Tam et al. (2013) Tam, Z. Y., J. Gruber, B. Halliwell, and R. Gunawan, 2013 Mathematical modeling of the role of mitochondrial fusion and fission in mitochondrial DNA maintenance. PloS ONE 8: e76230.
  * Tam et al. (2015) Tam, Z. Y., J. Gruber, B. Halliwell, and R. Gunawan, 2015 Context-dependent role of mitochondrial fusion-fission in clonal expansion of mtDNA mutations. PLoS Comp. Biol. 11: e1004183.
  * Taylor et al. (2003) Taylor, R. W., M. J. Barron, G. M. Borthwick, A. Gospel, P. F. Chinnery, et al., 2003 Mitochondrial DNA mutations in human colonic crypt stem cells. J. Clin. Invest. 112: 1351–1360.
  * Trifunovic et al. (2004) Trifunovic, A., A. Wredenberg, M. Falkenberg, J. N. Spelbrink, A. T. Rovio, et al., 2004 Premature ageing in mice expressing defective mitochondrial DNA polymerase. Nature 429: 417.
  * Twig et al. (2008) Twig, G., A. Elorza, A. J. Molina, H. Mohamed, J. D. Wikstrom, et al., 2008 Fission and selective fusion govern mitochondrial segregation and elimination by autophagy. EMBO J. 27: 433–446.
  * Van Kampen (1992) Van Kampen, N. G., 1992 Stochastic processes in physics and chemistry, volume 1. Elsevier.
  * Wachsmuth et al. (2016) Wachsmuth, M., A. Huebner, M. Li, B. Madea, and M. Stoneking, 2016 Age-related and heteroplasmy-related variation in human mtDNA copy number. PLoS Genet. 12: e1005939.
  * Wai et al. (2008) Wai, T., D. Teoli, and E. A. Shoubridge, 2008 The mitochondrial DNA genetic bottleneck results from replication of a subpopulation of genomes. Nat. Genet. 40: 1484.
  * Wallace and Chalkia (2013) Wallace, D. C. and D. Chalkia, 2013 Mitochondrial DNA genetics and the heteroplasmy conundrum in evolution and disease. Cold Spring Harb. Perspect. Biol. 5: a021220.
  * Wilkinson (2011) Wilkinson, D. J., 2011 Stochastic modelling for systems biology. CRC press.
  * Wonnapinij et al. (2008) Wonnapinij, P., P. Chinnery, and D. Samuels, 2008 The distribution of mitochondrial DNA heteroplasmy due to random genetic drift. Am. J. Hum. Genet. 83: 582.
  * Woods et al. (2018) Woods, D. C., K. Khrapko, and J. L. Tilly, 2018 Influence of maternal aging on mitochondrial heterogeneity, inheritance, and function in oocytes and preimplantation embryos. Genes 9: 265.
  * Ye et al. (2014) Ye, K., J. Lu, F. Ma, A. Keinan, and Z. Gu, 2014 Extensive pathogenicity of mitochondrial heteroplasmy in healthy human individuals. Proc. Natl. Acad. Sci. USA 111: 10654–10659.
  * Zamponi et al. (2018) Zamponi, N., E. Zamponi, S. A. Cannas, O. V. Billoni, P. R. Helguera, et al., 2018 Mitochondrial network complexity emerges from fission/fusion dynamics. Sci. Rep. 8: 363.
  * Zheng et al. (2006) Zheng, W., K. Khrapko, H. A. Coller, W. G. Thilly, and W. C. Copeland, 2006 Origins of human mitochondrial point mutations as DNA polymerase $\gamma$-mediated errors. Mutat. Res. 599: 11–20.

## Supporting Information

## 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

$\displaystyle s+s$ $\displaystyle\xrightarrow{\gamma}$ $\displaystyle f+f$
(S1) $\displaystyle s+f$ $\displaystyle\xrightarrow{\gamma}$ $\displaystyle
f+f$ (S2) $\displaystyle f$ $\displaystyle\xrightarrow{\beta}$ $\displaystyle
s$ (S3) $\displaystyle s$ $\displaystyle\xrightarrow{\alpha\rho}$
$\displaystyle s+s$ (S4) $\displaystyle f$
$\displaystyle\xrightarrow{\alpha\rho}$ $\displaystyle f+f$ (S5)
$\displaystyle s$ $\displaystyle\xrightarrow{\eta\rho}$
$\displaystyle\emptyset$ (S6) $\displaystyle f$
$\displaystyle\xrightarrow{\rho}$ $\displaystyle\emptyset$ (S7)

where Eq. (S1)-(S3) are analogous to Eq. (1)-(3) where mutant species are
neglected. Eq. (S4) and (S5) are simple birth processes with a shared constant
rate $\alpha\rho$. Eq. (S6) and (S7) are simple death processes with rates
$\eta\rho$ and $\rho$ respectively. The parameter $\rho$ is shared amongst all
of the birth and death reactions in Eqs. (S4)–(S7). $\rho$ 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 $\eta$. Since additional death is introduced into the system when
$\eta>1$, we include the parameter $\alpha>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

$\displaystyle\frac{\text{d}s}{\text{d}t}$ $\displaystyle=$
$\displaystyle-\gamma s^{2}-\gamma fs+\beta f+\alpha\rho s-\eta\rho s$ (S8)
$\displaystyle\frac{\text{d}f}{\text{d}t}$ $\displaystyle=$
$\displaystyle\gamma s^{2}+\gamma fs-\beta f+\alpha\rho f-\rho f$ (S9)

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
$\gamma$ rather than $2\gamma$ (justification of this is presented below, see
Eq. (S20)). 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 as 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

$s_{1j}X_{1}+\dots+s_{Nj}X_{N}\xrightarrow{\hat{k}_{j}}r_{1j}X_{1}+\dots+r_{Nj}X_{N}$
(S10)

where $s_{ij}$ and $r_{ij}$ are stoichiometric coefficients. We define
$\hat{k}_{j}$ as the microscopic rate for this reaction. The dimensionality of
this parameter will vary depending upon the stoichiometric coefficients
$s_{ij}$. $\hat{k}_{j}$ 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
$\mathbf{x}=(x_{1},\dots,x_{N})^{T}$, where $x_{i}$ is the copy number of the
$i$th species, allows us to write the CME as (Grima, 2010)

$\frac{\partial P(\mathbf{x},t|\mathbf{x}_{0},t_{0})}{\partial
t}=\Omega\sum_{j=1}^{R}\left(\prod_{i=1}^{N}E_{i}^{-S_{ij}}-1\right)\hat{f}_{j}(\mathbf{x},\Omega)P(\mathbf{x},t|\mathbf{x}_{0},t_{0})$
(S11)

where $\Omega$ 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_{i}^{-S_{ij}}$ is referred to as the step operator and is
defined through the relation
$E_{i}^{-S_{ij}}(g(\mathbf{x}))=g(x_{1},\dots,x_{i}-S_{ij},\dots,x_{N})$, for
any function of state $g(\mathbf{x})$. $\hat{f}_{j}(\mathbf{x},\Omega)$ is the
microscopic rate function of reaction $j$, which in general depends on both
the state and the system size. A factor of $\Omega$ 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
$\mathbf{x}$, the $j$th reaction occurs in the time interval $[t,t+\text{d}t)$
somewhere in $\Omega$ (Gillespie, 2007) is

$\hat{a}_{j}(\mathbf{x},\Omega)\text{d}t\vcentcolon=\Omega\hat{f}_{j}(\mathbf{x},\Omega)\text{d}t.$
(S12)

$\hat{a}_{j}(\mathbf{x},\Omega)$ is termed the propensity function (or
“hazard”) and is of particular relevance in the stochastic simulation
algorithm (Gillespie, 1976), since
$\hat{a}_{j}(\mathbf{x},\Omega)/\sum_{j}\hat{a}_{j}(\mathbf{x},\Omega)$
determines the probability that the $j$th reaction occurs next.

For the microscopic rate function, we may write

$\hat{f}_{j}(\mathbf{x},\Omega)=\hat{k}_{j}\prod_{i=1}^{N}\Omega^{-s_{ij}}\binom{x_{i}}{s_{ij}}.$
(S13)

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)

$\frac{\text{d}\phi_{i}}{\text{d}t}=\sum_{j=1}^{R}S_{ij}\tilde{f}_{j}(\boldsymbol{\phi})$
(S14)

where $\boldsymbol{\phi}=(\phi_{1},\dots,\phi_{N})^{T}$ is the vector of
macroscopic concentrations (of dimensions molecules per unit volume) and
$\tilde{f}_{j}(\boldsymbol{\phi})$ is the macroscopic rate function satisfying

$\tilde{f}_{j}(\boldsymbol{\phi})=\tilde{k}_{j}\prod_{i=1}^{N}\phi_{i}^{s_{ij}}$
(S15)

where $\tilde{k}_{j}$ is the macroscopic rate for the $j$th reaction. We
distinguish between $\hat{k}_{j}$ and $\tilde{k}_{j}$, 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 $\tilde{k}_{j}$ 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
Eq. (S11) to correspond to the macroscopic description in Eq. (S14), 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_{x_{i}\rightarrow\infty}\hat{f}_{j}(\mathbf{x},\Omega)=\tilde{f}_{j}(\boldsymbol{\phi})\quad\forall\
i,j.$ (S16)

In applying this constraint on all species $i$ and all reactions $j$, we may
derive a general relationship between $\hat{k}_{j}$ and $\tilde{k}_{j}$

$\lim_{x_{i}\rightarrow\infty}\hat{k}_{j}\prod_{i=1}^{N}\Omega^{-s_{ij}}\frac{x_{i}!}{s_{ij}!(x_{i}-s_{ij})!}=\tilde{k}_{j}\prod_{i=1}^{N}\phi_{i}^{s_{ij}}.$
(S17)

We can make two approximations to generate a more convenient relationship
between the microscopic and macroscopic rates. Firstly, we assume that

$x_{i}\approx\Omega\phi_{i}.$ (S18)

This is a small noise approximation, since it is often assumed that
$x_{i}=\Omega\phi_{i}+\Omega^{1/2}\xi_{i}$, where $\xi_{i}$ is a noise term
(Van Kampen, 1992). If $\xi_{i}$ is small then $x_{i}\approx\Omega\phi_{i}$ is
a valid approximation. Secondly, we assume that

$x_{i}(x_{i}-1)\dots(x_{i}-s_{ij}+1)\approx x_{i}^{s_{ij}}.$ (S19)

This is a large copy number approximation: in the case of e.g. a bimolecular
reaction ($2X_{i}\rightarrow*$) with $s_{ij}=2$, the approximation is of the
form $x_{i}(x_{i}-1)\approx x_{i}^{2}$ or $x_{i}\approx x_{i}-1$. By applying
Eq. (S19) to the factor of $x_{i}!$ in Eq. (S17), the factor of
$(x_{i}-s_{ij})!$ cancels from the left-hand side. Simplifying using Eq.
(S18), $\phi_{i}^{s_{ij}}$ cancels from both sides and we arrive at the
important relationship

$\tilde{k}_{j}\approx\hat{k}_{j}\prod_{i=1}^{N}\frac{1}{s_{ij}!}.$ (S20)

With Eq. (S14), Eq. (S15) and Eq. (S20) 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 $\hat{k}_{j}$. This equation highlights that
for reactions with $s_{ij}\geq 2$, $\tilde{k}_{j}\neq\hat{k}_{j}$, as is the
case for bimolecular reactions of the form $2X_{i}\rightarrow*$ (see Eq. (1)
and Eq. (S1)).

Importantly, if the microscopic rate function is a function of state then
$\hat{k}=\hat{k}(\mathbf{x})$ and
$\tilde{k}=\tilde{k}(\boldsymbol{\phi})\approx\tilde{k}(\mathbf{x}/\Omega)$.
In this case, Eq. (S20) still applies since the above argument assumed nothing
about the particular forms of $\hat{k}$ and $\tilde{k}$. However, additional
factors of $\Omega^{-1}$ are induced by applying Eq. (S18), which may carry
through to the individual parameters of $\tilde{k}(\boldsymbol{\phi})$. 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 Eq. (1)-Eq. (9) in the Main Text. In accordance with the notation
in the previous section, we will redefine all of the rates in Eq. (1)-Eq. (10)
with a hat notation ($\hat{a}$, for a general rate parameter $a$), to reflect
that these are stochastic rates. Deterministic rates will be denoted with a
tilde ($\tilde{a}$). Our aim will be to write a set of ODEs in terms of the
stochastic rates, $\hat{a}$, for which we are able to estimate values.

We will begin by considering the fusion network equations Eq. (1) and Eq. (2).
For clarity, we rewrite Eq. (1) to allow the reaction to proceed with some
arbitrary rate $\hat{\rho}$:

$X_{S}+X_{S}\xrightarrow{\hat{\rho}}X_{F}+X_{F},$ (S21)

where $X$ denotes either a wild-type ($W$) or mutant ($M$). We will
subsequently fix $\hat{\rho}$ to the rate of all other fusion reactions
$\hat{\gamma}$. We do this because Eq. (1) is a bimolecular reaction involving
one species: a fundamentally different reaction to bimolecular reactions
involving two species, as we will now see.

Since $\hat{\rho},\hat{\gamma}=$ const, we may use Eq. (S20), resulting in the
deterministic rates

$\displaystyle\tilde{\rho}=\frac{\hat{\rho}}{2}$ (S22)
$\displaystyle\tilde{\gamma}=\hat{\gamma}$ (S23)

for Eq. (S21) and Eq. (2) respectively. If we then enforce the microscopic
rates to be equal for both of these fusion reactions, i.e.
$\hat{\rho}=\hat{\gamma}$, then $\tilde{\rho}=\hat{\gamma}/2$. All other
fusion reactions have $\tilde{\gamma}=\hat{\gamma}$ by application of Eq.
(S20). Application of Eq. (S20) to the fission reaction in Eq. (3) shows that
$\hat{\beta}=\tilde{\beta}$.

For Eq. (4), we have chosen a $\hat{\lambda}$ which is not a constant, but a
function of the copy numbers of the chemical species
($\hat{\lambda}=\hat{\lambda}(\mathbf{x})$ where
$\mathbf{x}=(w_{s},w_{f},m_{s},m_{f})$, see Eq. (10)). As pointed out in the
previous section, care must be taken in writing down the deterministic
analogue of $\hat{\lambda}$. Applying Eq. (S20), we have

$\hat{\mu}+\hat{b}(\hat{\kappa}-(w_{s}+w_{f}+\hat{\delta}m_{s}+\hat{\delta}m_{f}))=\tilde{\mu}+\tilde{b}(\tilde{\kappa}-(\phi_{w_{s}}+\phi_{w_{f}}+\tilde{\delta}\phi_{m_{s}}+\tilde{\delta}\phi_{m_{f}})).$
(S24)

Applying Eq. (S18) and equating individual terms, we arrive at

$\displaystyle\hat{\mu}$ $\displaystyle=$ $\displaystyle\tilde{\mu}$ (S25)
$\displaystyle\hat{b}w_{s}=\tilde{b}\frac{w_{s}}{\Omega}$
$\displaystyle\implies$ $\displaystyle\hat{b}=\tilde{b}/\Omega$ (S26)
$\displaystyle\hat{b}\hat{k}=\tilde{b}\tilde{k}$ $\displaystyle\implies$
$\displaystyle\hat{k}=\tilde{k}\Omega$ (S27)
$\displaystyle\hat{b}\hat{\delta}m_{s}=\tilde{b}\tilde{\delta}\phi_{m_{s}}$
$\displaystyle\implies$ $\displaystyle\hat{\delta}=\tilde{\delta}.$ (S28)

In this study, we let $\Omega=1$ so the above 4 parameters are identical to
their deterministic counterparts. Hence, by application of Eq. (S14), we
arrive at the following set of ODEs

$\displaystyle\begin{split}\frac{\text{d}\phi_{w_{s}}}{\text{d}t}=-2\cdot\frac{\hat{\gamma}}{2}\phi_{w_{s}}^{2}-\hat{\gamma}\phi_{w_{s}}\\\
\phi_{w_{f}}+\hat{\beta}\phi_{w_{f}}-(\hat{\mu}+\hat{b}(\hat{\kappa}-(\phi_{w_{s}}+\phi_{w_{f}}\\\
+\hat{\delta}\phi_{m_{s}}+\hat{\delta}\phi_{m_{f}})))\phi_{w_{s}}-\hat{\mu}\phi_{w_{s}}-\\\
\hat{\gamma}\phi_{m_{f}}\phi_{w_{s}}-\hat{\gamma}\phi_{w_{s}}\phi_{m_{s}}\end{split}$
(S29)
$\displaystyle\begin{split}\frac{\text{d}\phi_{m_{s}}}{\text{d}t}=-2\cdot\frac{\hat{\gamma}}{2}\phi_{m_{s}}^{2}-\hat{\gamma}\phi_{m_{s}}\\\
\phi_{m_{f}}+\hat{\beta}\phi_{m_{f}}-(\hat{\mu}+\hat{b}(\hat{\kappa}-(\phi_{w_{s}}+\phi_{w_{f}}\\\
+\hat{\delta}\phi_{m_{s}}+\hat{\delta}\phi_{m_{f}})))\phi_{m_{s}}-\hat{\mu}\phi_{m_{s}}-\\\
\hat{\gamma}\phi_{w_{f}}\phi_{m_{s}}-\hat{\gamma}\phi_{w_{s}}\phi_{m_{s}}\end{split}$
(S30)
$\displaystyle\begin{split}\frac{\text{d}\phi_{w_{f}}}{\text{d}t}=2\cdot\frac{\hat{\gamma}}{2}\phi_{w_{s}}^{2}+\hat{\gamma}\phi_{w_{s}}\\\
\phi_{w_{f}}-\hat{\beta}\phi_{w_{f}}+(\hat{\mu}+\hat{b}(\hat{\kappa}-(\phi_{w_{s}}+\phi_{w_{f}}\\\
+\hat{\delta}\phi_{m_{s}}+\hat{\delta}\phi_{m_{f}})))(2\phi_{w_{s}}+\phi_{w_{f}})+\\\
\hat{\gamma}\phi_{m_{f}}\phi_{w_{s}}+\hat{\gamma}\phi_{w_{s}}\phi_{m_{s}}\end{split}$
(S31)
$\displaystyle\begin{split}\frac{\text{d}\phi_{m_{f}}}{\text{d}t}=2\cdot\frac{\hat{\gamma}}{2}\phi_{m_{s}}^{2}+\hat{\gamma}\phi_{m_{s}}\\\
\phi_{m_{f}}-\hat{\beta}\phi_{m_{f}}+(\hat{\mu}+\hat{b}(\hat{\kappa}-(\phi_{w_{s}}+\phi_{w_{f}}\\\
+\hat{\delta}\phi_{m_{s}}+\hat{\delta}\phi_{m_{f}})))(2\phi_{m_{s}}+\phi_{m_{f}})+\\\
\hat{\gamma}\phi_{w_{f}}\phi_{m_{s}}+\hat{\gamma}\phi_{w_{s}}\phi_{m_{s}}.\end{split}$
(S32)

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

$\displaystyle\begin{split}x_{1}=(b^{2}(\beta^{2}+2\beta(\gamma\kappa+3\mu)+\gamma^{2}\kappa^{2}+\mu^{2}+2\gamma\mu\\\
(\kappa+2(\delta-1)\phi_{m_{s}}))+2b\gamma\mu(-\beta+\mu+\gamma(\kappa-2\delta\\\
\phi_{m_{s}}+2\phi_{m_{s}}))+\gamma^{2}\mu^{2})^{1/2}\end{split}$ (S33)

the non-trivial, physically-realizable, component of the steady state may be
parametrized in terms of $\phi_{m_{s}}$ and written as

$\displaystyle\begin{split}\phi_{w_{s}}=-(\beta
b^{2}\kappa+b^{2}\gamma\kappa^{2}+b^{2}\kappa\mu+2b^{2}\delta\mu\\\
\phi_{m_{s}}+\beta^{2}b+\beta b\gamma\kappa+5\beta
b\mu+3b\gamma\kappa\mu-\beta\gamma\\\
\mu+2b\mu^{2}-2b\gamma\delta\mu\phi_{m_{s}}-2b\gamma\mu\phi_{m_{s}}-x_{1}(b\kappa+\\\
\beta+2\mu)+2\gamma\mu^{2}+2\gamma^{2}\mu\phi_{m_{s}})/(2\mu(b-\gamma)^{2})\end{split}$
(S34)
$\displaystyle\begin{split}\phi_{w_{f}}=(b^{3}(\beta^{2}\kappa+\beta(\gamma\kappa(\kappa+(\delta-1)\phi_{m_{s}})+\mu(3\\\
\kappa+\delta\phi_{m_{s}}))+\phi_{m_{s}}(\gamma^{2}(\delta-1)\kappa^{2}-\delta\mu^{2}+\gamma\mu\\\
((2\delta-3)\kappa+2(\delta-1)\delta\phi_{m_{s}})))+b^{2}(\beta^{3}+\beta^{2}(6\mu+\gamma\\\
(\kappa+(\delta-1)\phi_{m_{s}}))+\beta(5\mu^{2}+\gamma^{2}\delta\kappa\phi_{m_{s}}+\gamma^{2}\\\
\kappa(-\phi_{m_{s}})+6\gamma\delta\mu\phi_{m_{s}}-8\gamma\mu\phi_{m_{s}}-\kappa
x_{1})+\\\
\phi_{m_{s}}(2\gamma^{2}\mu(\delta\kappa+\delta^{2}(-\phi_{m_{s}})+\phi_{m_{s}})+\gamma((6\delta\\\
-5)\mu^{2}-(\delta-1)\kappa x_{1})-\delta\mu
x_{1}))-b(\beta(\gamma^{2}\mu(-(\kappa+(4-3\\\
\delta)\phi_{m_{s}}))+\gamma(2\mu^{2}+(\delta-1)\phi_{m_{s}}x_{1})+3\mu
x_{1})+\gamma\mu\\\
\phi_{m_{s}}(\gamma(\delta-2)\mu+\gamma^{2}(\kappa-2(\delta-1)\phi_{m_{s}})+(\delta-3)x_{1})+\\\
\
\beta^{2}(2\gamma\mu+x_{1}))+\gamma\mu(\beta-\gamma\phi_{m_{s}})(\gamma\mu+x_{1}))/(2b\\\
\mu(b-\gamma)^{2}(\beta+\gamma(\delta-1)\phi_{m_{s}}))\end{split}$ (S35)
$\displaystyle\begin{split}\phi_{m_{f}}=(\phi_{m_{s}}(-b\beta+b\gamma\kappa+b\mu-2b\gamma\delta\phi_{m_{s}}+2b\\\
\gamma\phi_{m_{s}}+\gamma\mu+x_{1}))/(2b(\beta+\gamma(\delta-1)\phi_{m_{s}})).\end{split}$
(S36)

Since the steady state is parametrized by $\phi_{m_{s}}$, the steady state is
therefore a line.

## Proof of heteroplasmy relation for linear feedback control

In this section we show that Eq. (13) holds for the system described by Eq.
(1)-Eq. (9) given the replication rate in Eq. (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 Eq. (S11) (CME) as

$\frac{\partial P(\mathbf{x},t)}{\partial
t}=\int_{-\infty}^{\infty}\operatorname{d}\\!{}\mathbf{x}{}^{\prime}\left[T(\mathbf{x}{}|\mathbf{x}{}^{\prime})P(\mathbf{x}{}^{\prime},t)-T(\mathbf{x}{}^{\prime}|\mathbf{x}{})P(\mathbf{x}{},t)\right]$
(S37)

where we have set $\Omega=1$. $T(\mathbf{x}{}|\mathbf{x}{}^{\prime})$ is the
transition rate from state $\mathbf{x}{}^{\prime}\rightarrow\mathbf{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

$\frac{\partial P(\mathbf{x},t)}{\partial
t}\approx\int_{-\infty}^{\infty}\left(-\nabla\left(T(\mathbf{x}^{\prime}|\mathbf{x})P(\mathbf{x})\right)^{T}\cdot(\mathbf{x}^{\prime}-\mathbf{x})+\frac{1}{2}(\mathbf{x}^{\prime}-\mathbf{x})^{T}\cdot\mathbf{H}\cdot(\mathbf{x}^{\prime}-\mathbf{x})\right)\operatorname{d}\\!{}\mathbf{x}^{\prime}$
(S38)

where $\mathbf{H}(\mathbf{x})$ is the Hessian matrix of
$T(\mathbf{x}^{\prime}|\mathbf{x})P(\mathbf{x})$

$\mathbf{H}\vcentcolon=\begin{pmatrix}\frac{\partial^{2}}{\partial
x_{1}^{2}}&\dots&\frac{\partial^{2}}{\partial x_{1}\partial x_{N}}\\\
\vdots&&\vdots\\\ \frac{\partial^{2}}{\partial x_{N}\partial
x_{1}}&\dots&\frac{\partial^{2}}{\partial
x_{N}^{2}}\end{pmatrix}T(\mathbf{x}^{\prime}|\mathbf{x})P(\mathbf{x})$ (S39)

(see (Gardiner, 1985) for a proof of this in the univariate case).

A transition to each possible neighbouring state $\mathbf{x}^{\prime}$
corresponds to some reaction $j$ which moves the state from
$\mathbf{x}\rightarrow\mathbf{x}^{\prime}$. Since we know the influence of
each reaction on state $\mathbf{x}$ through the constant stoichiometry matrix
$S_{ij}$, and that the propensity of a reaction does not depend upon
$\mathbf{x}^{\prime}$ itself (see Eq. (S13)), we may transition from a
notation involving $\mathbf{x}$ and $\mathbf{x}^{\prime}$ into a notation
involving $\mathbf{x}$ and $j$. We may therefore define
$T_{j}(\mathbf{x})\vcentcolon=T(\mathbf{x}^{\prime}|\mathbf{x})\equiv\hat{f}_{j}(\mathbf{x})$
(see Eq. (S13)), and let
$\mathbf{H}(\mathbf{x})\rightarrow\mathbf{H}_{j}(\mathbf{x})$.

We now make a large copy number assumption in order to simplify
$T_{j}(\mathbf{x}{})$. To take a large copy number limit, we assume that
$x_{i}!\approx x_{i}^{s_{ij}}(x_{i}-s_{ij})!$ resulting in

$T_{j}(\mathbf{x})\approx\hat{k}_{j}\prod_{i=1}^{N}\frac{x_{i}^{s_{ij}}}{s_{ij}!}.$
(S40)

This approximation is exact when $s_{ij}=0,1$, but inexact when $s_{ij}\geq
2$. For example, if we consider the second-order bimolecular reaction in Eq.
(1), Eq. (S40) is equivalent to assuming $w_{s}^{2}\approx w_{s}(w_{s}-1)$;
consequently, a factor of $1/(s_{ij}!)=1/2$ arises in $T_{j}(\mathbf{x})$ as a
combinatorial factor from stochastic considerations.

##### Fokker-Planck equation for chemical reaction networks

We now wish to re-write Eq. (S38) as a Fokker-Planck equation. Since the
integral in Eq. (S38) is over $\mathbf{x}^{\prime}$, and every
$\mathbf{x}^{\prime}$ corresponds to a reaction $j$, we may interpret the
integral in Eq. (S38) as a sum over all reactions, i.e.
$\int\operatorname{d}\\!{}\mathbf{x}^{\prime}\rightarrow\sum_{j=1}^{R}$.
Hence, for the $j$th reaction,
$[(\mathbf{x}^{\prime}-\mathbf{x})]_{i}=S_{ij}$. With these observations, we
may write the first integral of Eq. (S38) as

$\displaystyle\int_{-\infty}^{\infty}-\nabla(T(\mathbf{x}^{\prime}|\mathbf{x})P(\mathbf{x}))^{T}\cdot(\mathbf{x}^{\prime}-\mathbf{x})\operatorname{d}\\!{}\mathbf{x}^{\prime}$
$\displaystyle=\int_{-\infty}^{\infty}-\nabla(T_{j}(\mathbf{x})P(\mathbf{x},t))^{T}\cdot(\mathbf{x}^{\prime}-\mathbf{x})\operatorname{d}\\!{}\mathbf{x}^{\prime}$
$\displaystyle=-\sum_{j=1}^{R}\sum_{i=1}^{N}\frac{\partial}{\partial
x_{i}}(T_{j}(\mathbf{x})P(\mathbf{x},t))S_{ij}$
$\displaystyle=-\sum_{i=1}^{N}\frac{\partial}{\partial
x_{i}}A_{i}P(\mathbf{x},t)$ (S41)

where

$\mathbf{A}\vcentcolon=\mathbf{S}\cdot\mathbf{T}.$ (S42)

$\mathbf{A}$ is a vector of length $N$,
$[\mathbf{S}]_{ij}\vcentcolon=r_{ij}-s_{ij}$ is the $N\times R$ stoichiometry
matrix Eq. (S10), and $\mathbf{T}$ is the vector of transition rates, of
length $R$ (for which we have taken a large copy number approximation in Eq.
(S40)). To re-write the second integral of Eq. (S38), we write an element of
the Hessian $\mathbf{H}_{j}$ in Eq. (S39) as

$H_{jlm}=\frac{\partial^{2}}{\partial x_{l}\partial
x_{m}}T_{j}(\mathbf{x})P(\mathbf{x},t)$ (S43)

where $j=1,\dots,R$ and $l,m=1,\dots,N$. Thus, we may write

$\displaystyle\int_{-\infty}^{\infty}\frac{1}{2}(\mathbf{x}^{\prime}-\mathbf{x})^{T}\cdot\mathbf{H}_{j}\cdot(\mathbf{x}^{\prime}-\mathbf{x})\operatorname{d}\\!{}\mathbf{x}^{\prime}$
$\displaystyle=\frac{1}{2}\sum_{j=1}^{R}\sum_{l=1}^{N}\sum_{m=1}^{N}S_{lj}H_{jlm}S_{mj}$
$\displaystyle=\frac{1}{2}\sum_{j=1}^{R}\sum_{l=1}^{N}\sum_{m=1}^{N}S_{lj}\frac{\partial^{2}}{\partial
x_{l}\partial x_{m}}T_{j}P(\mathbf{x},t)S_{mj}$
$\displaystyle=\frac{1}{2}\sum_{l=1}^{N}\sum_{m=1}^{N}\frac{\partial^{2}}{\partial
x_{l}\partial
x_{m}}\left(\sum_{j=1}^{R}S_{lj}T_{j}S_{mj}\right)P(\mathbf{x},t)$
$\displaystyle=\frac{1}{2}\sum_{i,m=1}^{N}\frac{\partial^{2}}{\partial
x_{i}\partial x_{m}}B_{im}P(\mathbf{x},t)$ (S44)

where

$\mathbf{B}\vcentcolon=\mathbf{S}\cdot\operatorname{Diag}(\mathbf{T})\cdot\mathbf{S}^{T}.$
(S45)

$\mathbf{B}$ is an $N\times N$ matrix, and $\operatorname{Diag}(\mathbf{Y})$
is a diagonal matrix whose main diagonal is the vector $\mathbf{Y}$. We may
therefore re-write Eq. (S38) as a Fokker-Planck equation for the state vector
$\mathbf{x}$ of the form

$\frac{\partial P(\mathbf{x},t)}{\partial
t}\approx-\sum_{i=1}^{N}\frac{\partial}{\partial
x_{i}}[A_{i}(\mathbf{x})P(\mathbf{x},t)]+\frac{1}{2}\sum_{i,m=1}^{N}\frac{\partial^{2}}{\partial
x_{i}\partial x_{m}}[B_{im}(\mathbf{x})P(\mathbf{x})].$ (S46)

##### Fokker-Planck equation for an arbitrary function of state

We now wish to make a change of variables in Eq. (S46) to write down a Fokker-
Planck equation for an arbitrary scalar function of state $\mathbf{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 Eq. (S46) is equivalent (Jacobs, 2010) to the
following Itô stochastic differential equation (SDE)

$\operatorname{d}\\!{}\mathbf{x}=\mathbf{A}\operatorname{d}\\!{}t+\mathbf{G}\operatorname{d}\\!{}\mathbf{W}$
(S47)

where $\mathbf{G}\mathbf{G}^{T}\equiv\mathbf{B}$ (where $\mathbf{G}$ is an
$N\times R$ matrix) and $\operatorname{d}\\!{}\mathbf{W}$ is a vector of
independent Wiener increments of length $R$, and a Wiener increment
$\operatorname{d}\\!{W}$ satisfies

$\int_{0}^{t}\operatorname{d}\\!{W}\vcentcolon=W(t),\
P(W,t)\equiv\frac{1}{\sqrt{2\pi t}}e^{-W^{2}/(2t)}.$ (S48)

Itô’s formula states that, for an arbitrary function $h(\mathbf{x},t)$ where
$\mathbf{x}$ satisfies Eq. (S47), we may write the following SDE

$\operatorname{d}\\!{h}(\mathbf{x},t)=\left\\{\frac{\partial h}{\partial
t}+\left(\nabla
h\right)^{T}\mathbf{A}+\frac{1}{2}\operatorname{Tr}\left[\mathbf{G}^{T}\mathbf{H}_{h}(\mathbf{x})\mathbf{G}\right]\right\\}\operatorname{d}\\!{t}+(\nabla
h)^{T}\mathbf{G}\operatorname{d}\\!{}\mathbf{W},$ (S49)

where $\mathbf{H}_{h}(\mathbf{x})$ is the Hessian matrix of $h(\mathbf{x},t)$
(see Eq. (S39), where $T(\mathbf{x}^{\prime}|\mathbf{x})P(\mathbf{x})$ should
be replaced with $h(\mathbf{x},t)$). Given the form of $\mathbf{B}$ in Eq.
(S45) we let

$\mathbf{G}=\mathbf{S}\cdot\operatorname{Diag}(\sqrt{\mathbf{T}}),$ (S50)

which satisfies $\mathbf{G}\mathbf{G}^{T}\equiv\mathbf{B}$.

For convenience, we may also perform the transformation purely at the level of
Fokker-Planck equations. Let $h(\mathbf{x},t)$ satisfy the general Fokker-
Planck equation

$\frac{\partial P(h,t)}{\partial t}=-\frac{\partial}{\partial
h}[\tilde{A}(h,t)P(h,t)]+\frac{1}{2}\frac{\partial^{2}}{\partial
h^{2}}[\tilde{B}(h,t)P(h,t)]$ (S51)

for scalar functions $\tilde{A}(h,t)$ and $\tilde{B}(h,t)$. Using the cyclic
property of the trace in Eq. (S49), we may identify

$\tilde{A}=\frac{\partial h}{\partial t}+\left(\nabla
h\right)^{T}\mathbf{A}+\frac{1}{2}\operatorname{Tr}\left[\mathbf{B}H_{h}(\mathbf{x})\right]$
(S52)

where $\operatorname{Tr}$ is the trace operator. Also, from Eq. (S49),

$\tilde{B}=[(\nabla h)^{T}\mathbf{G}][(\nabla h)^{T}\mathbf{G}]^{T}=(\nabla
h)^{T}\mathbf{B}(\nabla h).$ (S53)

Hence, using Eq. (S51), Eq. (S52) and Eq. (S53), we may write down a Fokker-
Planck equation for an arbitrary function of state in terms of $\mathbf{A}$
and $\mathbf{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 Eq. (13). In what follows, we will
assume that $\mathbf{x}(t=0)=\mathbf{x}_{\text{ss}}$, where
$\mathbf{x}_{\text{ss}}$ is the state which is the solution of
$\mathbf{A}=\mathbf{0}$ (which is equivalent to finding the steady state
solution of the deterministic rate equation in Eq. (S14) due to our assumption
of large copy numbers and $\Omega=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 Eqs. (S34)–(S36)).
Inspection of the steady state solution reveals that the steady state depends
on the fusion ($\gamma$) and fission ($\beta$) 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 $\mathbb{V}(h)$.
We redefine the fusion and fission rates such that

$\displaystyle\gamma$ $\displaystyle\rightarrow$ $\displaystyle M\gamma$
$\displaystyle\beta$ $\displaystyle\rightarrow$ $\displaystyle M\beta$ (S54)

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

$h(\mathbf{x},t)=h(\mathbf{x})\vcentcolon=(m_{s}+m_{f})/(w_{s}+w_{f}+m_{s}+m_{f}),$
(S55)

as our choice for the function of state in the Fokker-Planck equation in Eq.
(S51). We will first compute the diffusion term $\tilde{B}$ for heteroplasmy
using Eq. (S53). If we constrain the state $\mathbf{x}$ to be forced onto the
steady state line $\mathbf{x}_{\text{ss}}$ (as per (Constable et al., 2016;
Parsons and Rogers, 2017)) in the high-churn limit, then upon defining

$\theta\vcentcolon=\sqrt{b^{2}(\beta+\gamma\kappa)^{2}-2b\gamma\mu\left(\beta-\gamma\kappa+2\gamma(\delta-1)m_{s}\right)+\gamma^{2}\mu^{2}}$
(S56)

we have

$\displaystyle\begin{split}\lim_{M\rightarrow\infty}\left(\left.\tilde{B}\right|_{\mathbf{x}=\mathbf{x}_{\text{ss}}}\right)=(16b^{2}\gamma^{2}\mu\\\
m_{s}(\beta+\gamma(\delta-1)m_{s}){}^{2}(b^{2}(\beta^{3}+2\beta^{2}\gamma\\\
(\delta-1)m_{s}+\beta\gamma^{2}(\kappa^{2}+(1-2\delta)\kappa
m_{s}+(\delta-1)(2\delta-1)\\\ m_{s}^{2})+\gamma^{3}\kappa
m_{s}((\delta-1)m_{s}-\kappa))+b(-\beta^{2}\\\
(2\gamma\mu+\theta)+\beta\gamma(\kappa(2\gamma\mu+\theta)+m_{s}(\gamma(5-4\\\
\delta)\mu-2(\delta-1)\theta))+\gamma^{2}m_{s}((\delta-1)m_{s}(3\gamma\mu\\\
+\theta)-\kappa(2\gamma\mu+\theta)))+\gamma\mu(\gamma\mu+\theta)\\\
(\beta-\gamma m_{s})))/(\beta^{3}(b(\gamma(\kappa-2\delta m_{s}\\\
+2m_{s})-\beta)+\gamma\mu+\theta){}^{4}).\end{split}$ (S57)

Eq. (S57) is difficult to understand. In order to perform further
simplification, we make an ansatz for the form of $\tilde{B}$
($\tilde{B}_{\text{An}}$) and seek to determine whether our ansatz is
equivalent to the derived form of $\tilde{B}$ under the constraints defined on
the left-hand side of Eq. (S57). Our ansatz takes the form

$\tilde{B}_{\text{An}}\vcentcolon=\lim_{M\rightarrow\infty}\left(\left.\frac{2\mu
h(1-h)}{n(\mathbf{x})}\cdot
f_{s}(\mathbf{x})\right|_{\mathbf{x}=\mathbf{x}_{\text{ss}}}\right)$ (S58)

where $f_{s}(\mathbf{x})\vcentcolon=(w_{s}+m_{s})/(w_{s}+w_{f}+m_{s}+m_{f})$
and $n(\mathbf{x})\vcentcolon=w_{s}+w_{f}+m_{s}+m_{f}$. Notice that this
ansatz is more general than Eq. (S57), since it has no explicit dependence
upon the parameters of the control law assumed in Eq. (10), and only
explicitly depends upon functions of state $\mathbf{x}$.

Upon substituting the steady state $\mathbf{x}_{\text{ss}}$ into the ansatz in
Eq. (S58) and taking the high-churn limit, we find that

$\displaystyle\begin{split}\tilde{B}_{\text{An}}=-(m_{s}(\beta+\gamma(\delta-1)m_{s})\\\
(b(\beta+\gamma\kappa)-\gamma\mu+\theta)\\\
(b(\beta+\gamma\kappa)+\gamma\mu-\theta)(b(m_{s}(2\beta\delta-\beta\\\
+\gamma\kappa)-2\beta\kappa)+m_{s}(\theta-\gamma\mu)))/(2b\beta^{3}(\kappa-\delta\\\
m_{s}+m_{s})^{2}(b(\gamma(\kappa-2\delta
m_{s}+2m_{s})-\beta)+\gamma\mu+\theta)).\end{split}$ (S59)

After some algebra (see GitHub repository for Mathematica notebook), it can be
shown that Eq. (S57) and Eq. (S59) are equivalent, i.e.

$\tilde{B}_{\text{An}}\equiv\lim_{M\rightarrow\infty}\left(\left.\tilde{B}\right|_{\mathbf{x}=\mathbf{x}_{\text{ss}}}\right).$
(S60)

As such, we may use $\tilde{B}$ and $\tilde{B}_{\text{An}}$ 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.

$\left.\tilde{A}\right|_{\mathbf{x}=\mathbf{x}_{\text{ss}}}\equiv 0.$ (S61)

A similar result is shown in (Constable et al., 2016) ((Equation S59)
therein). Substituting $h(\mathbf{x},t)=h$, $\tilde{A}$ and $\tilde{B}$ into
the Fokker-Planck equation for an arbitrary function of state Eq. (S51), we
have

$\frac{\partial P(h,t)}{\partial t}=\frac{1}{2}\frac{\partial^{2}}{\partial
h^{2}}\left[\left.\left(\frac{2\mu h(1-h)}{n(\mathbf{x})}\cdot
f_{s}(\mathbf{x})\right)\right|_{\mathbf{x}=\mathbf{x}_{\text{ss}}(h)}P(h,t)\right]$
(S62)

which is equivalent to the following SDE for heteroplasmy

$\operatorname{d}\\!{h}=\left.\sqrt{\frac{2\mu
h(1-h)f_{s}(\mathbf{x})}{n(\mathbf{x})}}\right|_{\mathbf{x}=\mathbf{x}_{\text{ss}}(h)}\operatorname{d}\\!{W}$
(S63)

in the limit of large network churn, large copy numbers, and a second-order
truncation of the Kramers-Moyal expansion. 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 $\mathbf{x}_{\text{ss}}(h)$ for
every displacement due to Wiener noise in $h$, and substitute into
$f_{s}(\mathbf{x})$ and $n(\mathbf{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 Eq. (S63)
may be approximated as constant. Since the general solution of the SDE

$\operatorname{d}\\!{y}=\sqrt{B}\operatorname{d}\\!{W}$ (S64)

for $B=$ const is

$y\sim\mathcal{N}(y|y_{0},Bt)$ (S65)

where $\mathcal{N}(y|y_{0},\sigma^{2})$ is a Gaussian distribution on $y$ with
mean $y_{0}$ and variance $\sigma^{2}$, and $y_{0}=y(t=0)$. Since we have
assumed that the state is initialised at
$\mathbf{x}(t=0)=\mathbf{x}_{\text{ss}}$, there are no deterministic transient
dynamics, so we may write

$\mathbb{V}(h)\approx\left.\frac{2\mu
t}{n(\mathbf{x})}h(\mathbf{x})(1-h(\mathbf{x}))f_{s}(\mathbf{x})\right|_{\mathbf{x}=\mathbf{x}_{\text{ss}}},$
(S66)

where $\mathbb{V}$ returns variance of a random variable. In this equation, we
take $\mathbf{x}=\mathbf{x}_{\text{ss}}=$ const, since we have assumed a low-
diffusion limit. We observe 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 Eq.
(S63). Firstly, we may consider replacing Eq. (4) with

$X_{S}\xrightarrow{\lambda}X_{S}+X_{S}.$ (S67)

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 Eq. (S63), 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

$X_{F}\xrightarrow{\xi\mu}\emptyset$ (S68)

where $0\leq\xi\leq 1$. We were not able to make analogous analytical progress
in this instance. However, numerical investigation (Figure S3E) revealed that
the following ansatz was able to predict heteroplasmy variance dynamics

$\mathbb{V}(h)\approx\left.\frac{2\mu
t}{n(\mathbf{x})}h(\mathbf{x})(1-h(\mathbf{x}))(f_{s}(\mathbf{x})+\xi(1-f_{s}(\mathbf{x})))\right|_{\mathbf{x}=\mathbf{x}_{\text{ss}}}.$
(S69)

In other words, allowing degradation of fused species results in a linear
correction to our heteroplasmy variance formula in Eq. (13). If fused species
are susceptible to degradation at the same rate as unfused species ($\xi=1$),
then $\mathbb{V}(h)$ loses $f_{s}$ dependence entirely and the mitochondrial
network has no influence over heteroplasmy dynamics.

We also explored various different forms of $\lambda(\mathbf{x})$ and
$\mu(\mathbf{x})$, which we label A-G after (Johnston and Jones, 2016), and
X-Z for several newly-considered functional forms, see Table S3 and Figure
S4A-I. Control D of (Johnston and Jones, 2016) involves no feedback, which we
do not explore – see Figure S1, and the discussion surrounding Eq. (S1). 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
$\mathbf{x}_{\text{ss}}$ in Eq. (S57), Eq. (S59) and Eq. (S61). For controls
B, C, E, F, G, Y and Z in Table S3, 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 Eq. (S63); this can be
shown numerically for controls A and X. However, controls G, Y and Z satisfy

$\operatorname{d}\\!{h}=\left.\sqrt{\frac{2\mu
h(1-h)}{n(\mathbf{x})}}\right|_{\mathbf{x}=\mathbf{x}_{\text{ss}}(h)}\operatorname{d}\\!{W}.$
(S70)

Notably, Eq. (S70) does not depend on $f_{s}$, unlike Eq. (S63) (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 Eq. (1)-Eq. (9), given the replication rate in Eq. (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$\pm$0.003 fissions/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 $\beta=33.12$
day-1.

The dimensions of $\beta$ are day-1 and not mitochondrion-1 day-1. This is
because if the propensity (see Eq. (S12), where $\Omega=1$) of e.g. Eq. (3) is
$\hat{a}_{\text{fis},w}=\beta w_{f}$ then the mean time to the next event is
$1/(\beta w_{f})$; therefore the dimension of $\beta$ 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 Eq. (S20)), each fusion reaction proceeds at rate
$\hat{a}_{\text{fus},j}\approx\gamma x_{i}^{2}$. Since there are 7 fusion
reactions (Eq. (1), Eq. (2), Eq. (7)-Eq. (9)), the total fusion propensity is
$\hat{a}_{\text{fus}}\approx 7\gamma x_{i}^{2}$. Similarly, the total fission
propensity is $\hat{a}_{\text{fis}}=\beta(w_{f}+m_{f})=2\beta 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,
$\hat{a}_{\text{fus}}=\hat{a}_{\text{fis}}$, and rearrange for the fusion rate
$\gamma$ to yield $\gamma=2\beta x_{i}/(7x_{i}^{2})\approx 3.8\times 10^{-2}$
day-1. The orders of magnitude difference between $\beta$ and $\gamma$ 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 ($\mu$) and the parameters of
the linear feedback control ($\kappa$, $b$ and $\delta$, see Eq. (10)). mtDNA
half-life is observed to be highly variable: in mice this can be between
10-100 days (Burgstaller et al., 2014a). 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 $\delta$ in the replication feedback control (see Eq. (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
$\delta=1$ as the simplest case of a neutral mutation (although $\delta\neq 1$
still defines a neutral model, since both mutant and wild-type alleles
experience the same replication and degradation rates per molecule, see
Eqs.(4)–(6)).

We are finally left with setting the parameters $\kappa$ and $b$ in the linear
feedback control Eq. (10). In the absence of a network state and mutants,
$\kappa$ is precisely equal to the steady state copy number, since the
degradation rate equals the replication rate when $\kappa=w$. However, the
presence of a network means that a subpopulation of mtDNAs (namely the fused
species) are immune to death, resulting in $\kappa$ 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 $\kappa$ 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 $\kappa$
could not be found whose deterministic steady state was sufficiently close to
$n=1000$. This result is intuitive because in the limit of $b\rightarrow 0$,
$\lambda=$ const. From the analysis above we have shown that constant genetic
rates ($\mu,\lambda$) 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\approx b^{*}$, and the corresponding
$\kappa$, such that the steady state copy number is controlled as weakly as
possible given the model structure.

## Rate renormalization

In Eqs. (1)–(10) we have neglected reactions such as

$X_{F}+X_{F}\rightarrow X_{F}+X_{F}$ (S71)

because they do not change the number of molecules in our state vector
$\mathbf{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. Eq. (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 Eq. (6)).
Therefore, if fission reactions always result in mitochondria containing