Mitochondrial network state scales mtDNA genetic dynamics
Juvid Aryaman, Charlotte Bowles, Nick S. Jones, Iain G. Johnston
MMitochondrial network state scales mtDNA genetic dynamics
Juvid Aryaman, ∗ , † , ‡ Charlotte Bowles, § Nick S. Jones, ∗ , ∗∗ , Iain G. Johnston, †† , §§ , ∗ 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 [email protected] [email protected] Abstract
Mitochondrial DNA (mtDNA) mutations cause severe congenital diseases but may also be associated withhealthy aging. MtDNA is stochastically replicated and degraded, and exists within organelles which undergodynamic fusion and fission. The role of the resulting mitochondrial networks in the time evolution of thecellular proportion of mutated mtDNA molecules (heteroplasmy), and cell-to-cell variability in heteroplasmy(heteroplasmy variance), remains incompletely understood. Heteroplasmy variance is particularly importantsince it modulates the number of pathological cells in a tissue. Here, we provide the first wide-reachingtheoretical framework which bridges mitochondrial network and genetic states. We show that, under a rangeof conditions, the (genetic) rate of increase in heteroplasmy variance and de novo mutation are proportionallymodulated 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 clearanceof mtDNA mutants. Our findings imply that modulating network state, mitophagy rate and copy number toslow down heteroplasmy dynamics when mean heteroplasmy is low could have therapeutic advantages formitochondrial disease and healthy aging.
Introduction
Mitochondrial DNA (mtDNA) encodes elements of the respiratory system vital for cellular function. Mutationof mtDNA is one of several leading hypotheses for the cause of normal aging (Kauppila et al. , 2017; L´opez-Ot´ın et al. , 2013), 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 encodescrucial 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 populationof mtDNA molecules (Aryaman et al. , 2019; Johnston, 2018; Stewart and Chinnery, 2015; Wallace and Chalkia,2013). The broad biomedical implications of mitochondrial DNA mutation, combined with the countablenature of mtDNAs and the stochastic nature of their dynamics, offer the opportunity for mathematicalunderstanding 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 oftentolerate relatively high levels of mtDNA mutation, until the fraction of mutated mtDNAs (termed heteroplasmy)exceeds a certain critical value where a pathological phenotype occurs (Aryaman et al. , 2017; Picard et al. ,2014; Rossignol et al. , 2003; Stewart and Chinnery, 2015). Fluctuations within individual cells mean that thefraction of mutant mtDNAs per cell is not constant within a tissue (Figure 1A), but follows a probabilitydistribution which changes with time (Figure 1B). Here, motivated by a general picture of aging, we will1 a r X i v : . [ q - b i o . S C ] J u l argely focus on the setting of non-dividing cells, which possess two mtDNA variants (although we willalso consider de novo mutation using simple statistical genetics models). The variance of the distributionof heteroplasmies gives the fraction of cells above a given pathological threshold (Figure 1B). Thereforeheteroplasmy variance is related to the number of dysfunctional cells above a phenotypic threshold withina tissue, and both heteroplasmy mean and variance are directly related to tissue physiology. Increases inheteroplasmy variance also increase the number of cells below a given threshold heteroplasmy, which can beadvantageous in e.g. selecting low-heteroplasmy embryos in pre-implantation genetic diagnosis for treatingmitochondrial disease (Burgstaller et al. , 2014b; Johnston et al. , 2015).Mitochondria exist within a network which dynamically fuses and fragments. Although the function ofmitochondrial networks remains an open question (Hoitzing et al. , 2015), it is often thought that a combinationof network dynamics and mitochondrial autophagy (termed mitophagy) act in concert to perform qualitycontrol on the mitochondrial population (Aryaman et al. , 2019; Johnston, 2018; Twig et al. , 2008). Observationsof 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 suggestedthat certain mtDNA mutations, such as deletions (Kowald and Kirkwood, 2018, 2013, 2014) and somepoint mutations (Li et al. , 2015; Lieber et al. , 2019; Samuels et al. , 2013; Ye et al. , 2014), are under theinfluence of selective effects. However, genetic models without selection have proven valuable in explainingthe 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 geneticallydistant (Burgstaller et al. , 2014a; Pan et al. , 2019), suggesting that neutral genetic theory may be usefulin 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; Johnston and Jones, 2016; Poovathingal et al. , 2009), potentially allowing us tobetter understand and quantify when selection is present. There is thus a set of open questions about howthe physical dynamics of mitochondria affect the genetic populations of mtDNA within and between cellsunder neutral dynamics.A number of studies have attempted to understand the impact of the mitochondrial network on mitochon-drial dysfunction through computer simulation (reviewed in Kowald and Klipp (2014)). These studies havesuggested: that clearance of damaged mtDNA can be assisted by high and funcitonally-selective mitochondrialfusion, or by intermediate fusion and selective mitophagy (Mouli et al. , 2009); that physical transport ofmitochondria 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 mitochondrialquality (Figge et al. , 2012).Despite providing valuable insights, these previous attempts to link mitochondrial genetics and networkdynamics, while important for breaking ground, have centered around complex computer simulations, makingit difficult to deduce general laws and principles. Here, we address this lack of a general theoretical frameworklinking 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 intuitiveunderstanding, 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 modellinghave shown success in understanding mitochondrial physiology from a purely genetic perspective (Capps et al. ,2003; Chinnery and Samuels, 1999; Johnston and Jones, 2016). Furthermore, there currently exists limitedevidence for pronounced, universal, selective differences of mitochondrial variants in vivo (Hoitzing, 2017;Stewart and Larsson, 2014). 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 themitochondrial network (no selective fusion). Evidence for negative selection of particular mtDNA mutationshas been observed in vivo (Morris et al. , 2017; Ye et al. , 2014); we therefore extend our analysis to exploreselectivity 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 mitochondrialgenetic and network dynamics, yet we find that a simple link emerges. For a broad range of situations, theexpansion of mtDNA mutants is strongly modulated by network state, such that the rate of increase of2 igure 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 atissue exhibit inter-cellular heterogeneity in mitochondrial content: both mutant load (heteroplasmy) and copynumber. ( B ) Inter-cellular heterogeneity implies that heteroplasmy is described by a probability distribution. Cellsabove a threshold heteroplasmy ( h ∗ , black dashed line) are thought to exhibit a pathological phenotype. Thelow-variance distribution (black line) has fewer cells above a pathological threshold heteroplasmy than thehigh-variance distribution (red line). Heteroplasmy is depicted as an approximately normal distribution, as this is theregime in which our approximations below hold: i.e. when the probability of fixation is small. ( C ) The chemicalreaction network we use to model the dynamics of mitochondrial DNA (see Main Text for a detailed description).MtDNAs are assigned a genetic state: mutant ( M ) or wild-type ( W ), and a network state: singleton (i.e. unfused, S )or fused ( F ). ( D ) The central result of our work is, assuming that a cell at time t = 0 is at its (deterministic)steady-state, heteroplasmy variance ( V ( h )) approximately increases with time ( t ), mitophagy rate ( µ ) and the fractionof mitochondria that are unfused ( f s ), and decreases with mtDNA copy number ( n ). Importantly, V ( h ) does notdepend on the absolute magnitude of the fission-fusion rates. Also see Table S1 for a summary of our key findings. de novo mutation, is proportional to the fractionof unfused mitochondria. We discover that this result stems from the general notion that fusion shieldsmtDNAs from turnover, since autophagy of large fragments of the mitochondrial network are unlikely, andconsequently rescales time. Importantly, we used our model for network dynamics to show that heteroplasmyvariance is independent of the absolute magnitude of the fusion and fission rates due to a separation oftimescales between genetic and network processes (in contrast to Tam et al. (2015)). Surprisingly, we findthe dependence of heteroplasmy statistics upon network state arises when the mitochondrial population sizeis controlled through replication, and vanishes when it is controlled through mitophagy, shedding new lighton 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 effectiveelimination of mitochondrial mutants (in contrast to Mouli et al. (2009)). We also confirm that mitophagyand mitochondrial DNA copy number affect the rate of accumulation of de novo mutations (Johnston andJones, 2016), see Table S1 for a summary of our key findings. We suggest that pharmacological interventionswhich promote fusion, slow mitophagy and increase copy number earlier in development may slow the rate ofaccumulation of pathologically mutated cells, with implications for mitochondrial disease and aging. Materials and Methods
Stochastic modelling of the coupling between genetic and network dynamics ofmtDNA populations
Our modelling approach takes a chemical master equation perspective by combining a general model of neutralgenetic drift (for instance, see Chinnery and Samuels (1999); Johnston and Jones (2016)) with a model ofmitochondrial network dynamics. We seek to understand the influence of the mitochondrial network uponmitochondrial genetics. The network state itself is influenced by several factors including metabolic poise andthe respiratory state of mitochondria (Hoitzing et al. , 2015; Mishra and Chan, 2016; Szabadkai et al. , 2006),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 undergoingturnover (or “relaxed replication” Stewart and Chinnery (2015)). MtDNAs exist within mitochondria, whichundergo 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 themitochondrial network in a simple way, without the need to resort to a spatial model or consider the precisenetwork structure, allowing us to make analytic progress and derive interpretable formulae in a more generalrange of situations.Our model can be decomposed into three notional blocks (Figure 1C). Firstly, the principal networkprocesses denote fusion and fission of mitochondria containing mtDNAs of the same allele X S + X S γ −→ X F + X F (1) X F + X S γ −→ X F + X F (2) X F β −→ X S (3)where X denotes either a wild-type ( W ) or a mutant ( M ) mtDNA (therefore a set of chemical reactionsanalogous to Eq. (1)-(3) exist for both DNA species). γ and β are the stochastic rate constants for fusion andfission respectively.Secondly, mtDNAs are replicated and degraded through a set of reactions termed genetic processes. Acentral assumption is that all degradation of mtDNAs occur through mitophagy, and that only small pieces ofthe mitochondrial network are susceptible to mitophagy; for parsimony we take the limit of only the singletonsbeing susceptible to mitophagy X S λ −→ X F + X F (4) X F λ −→ X F + X F (5) X S µ −→ ∅ (6)where λ and µ are the replication and mitophagy rates respectively, which are shared by both W and M resulting in a so-called ‘neutral’ genetic model. Eq. (6) denotes removal of the species from the system. The4ffect 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 singletonswhich 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 beconsidered neutral since both W and M possess the same replication and degradation rates per molecule ofmtDNA at any instance in time.Finally, mtDNAs of different genotypes may interact through fusion via a set of reactions we term networkcross-processes: W F + M S γ −→ W F + M F (7) M F + W S γ −→ M F + W F (8) W S + M S γ −→ W F + M F . (9)Any fusion or fission event which does not involve the generation or removal of a singleton leaves our systemunchanged; we term such events as non-identity-changing processes, which can be ignored in our system (seeSupporting Information, Rate renormalization for a discussion of rate renormalization). We have neglected denovo mutation in the model description above (although we will consider de novo mutation using a modifiedinfinite sites Moran model below).We found that treating λ = 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 iscontrolled to a particular value, as has been done by previous authors (Capps et al. , 2003; Chinnery andSamuels, 1999; Johnston and Jones, 2016). Allowing lower-case variables to denote the copy number of theirrespective molecular species, we will focus on a linear replication rate of the form (Hoitzing, 2017; Hoitzing et al. , 2017): λ = λ ( w T , m T ) = µ + b ( κ − ( w T + δ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 copynumber, and κ is a parameter which is indicative of (but not equivalent to) the steady state copy number. δ indicates the relative contribution of mutant mtDNAs to the control strength and is linked to the “maintenanceof wild-type” hypothesis (Durham et al. , 2007; Stewart and Chinnery, 2015). When 0 ≤ δ <
1, and bothmutant 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 mutantshave 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. λ may be re-written as λ = k + k w T + k m T for constants k i , and so only consists of 3 independent parameters. However we will retain λ in the formof Eq. (10) since the parameters µ , b , κ , and δ have the distinct physiological meanings described above(Hoitzing, 2017; Hoitzing et al. , 2017). Furthermore, λ may in general also depend on other cellular featuressuch as mitochondrial reactive oxygen species. Here, we seek to explain mitochondrial behaviour under asimple set of governing principles, but our approach can naturally be combined with a description of theseadditional factors to build a more comprehensive model. Analogues of this model (without a network) havebeen applied to mitochondrial systems (Capps et al. , 2003; Chinnery and Samuels, 1999). Overall, our simplemodel consists of 4 species ( W S , W F , M S , M F ), 6 independent parameters and 15 reactions, and captures thecentral 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 (Aryaman et al. , 2019; Samuels et al. , 2010; Wonnapinij et al. , 2008): h ( x ) = ( m s + m f ) / ( w s + w f + m s + m f ) (11)where x = ( w s , w f , m s , m f ) is the state of the system (not to be confused with mitochondrial “respiratorystates”). Hence, a heteroplasmy of h = 1 denotes a cell with 100% mutant mtDNA (i.e. a homoplasmic cell inthe mutant allele). Arguably, “mutant allele fraction” would be a more precise description of Eq.(11) but weretain the use of heteroplasmy for consistency. To convert to a definition of heteroplasmy which is maximalwhen the mutant allele fraction is 50%, one may simply use the conversion 0 . − | h ( x ) − . | .5 tatistical Analysis In Figures S3B, S4A-I, we compare Eq. (13) and Eq. (S72) to stochastic simulations, for various parametriza-tions and replication/degradation rates. To quantify the accuracy of these equations in predicting V ( h, t ), wedefine the following error metric (cid:15) (cid:15) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ˙ V ( h, t ) Th E t ( ˙ V ( h, t ) Sim ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (12)where ˙ V ( h, t ) is the time derivative of heteroplasmy variance with subscripts denoting theory (Th) andsimulation (Sim). An expectation over time ( E t ) is taken for the stochastic simulations, whereas ˙ V ( h, t ) is ascalar 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 varianceover time, independently of fission-fusion rate magnitudes
We first performed a deterministic analysis of the system presented in Eqs. (1)–(10), by converting the reactionsinto an analogous set of four coupled ordinary differential equations (see Eqs. (S29)–(S32)), and choosinga 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 thatcopy 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 reachingthis line, total copy number remains constant (Figure S2A) and the state of the system ceases to change withtime. This is a consequence of performing a deterministic analysis, which neglects stochastic effects, andour choice of replication rate in Eq. (10) which decreases with total copy number when w T + δm T > κ andvice versa, guiding the total population to a fixed total copy number. Varying the fission ( β ) and fusion ( γ )rates revealed a negative linear relationship between the steady-state fraction of singletons and copy number(Figure S2B).We may also simulate the system in Eqs. (1)–(9) stochastically, using the stochastic simulation algorithm(Gillespie, 1976), which showed that mean copy number is slightly perturbed from the deterministic predictiondue to the influence of variance upon the mean (Grima et al. , 2011; Hoitzing, 2017) (Figure 2C). Thestationarity of total copy number is a consequence of using δ = 1 for our nominal parametrization (i.e. theline of steady states is also a line of constant copy number). Choosing δ (cid:54) = 1 results in a difference in carryingcapacities between the two species, and non-stationarity of mean total copy number, as trajectories spreadalong the line of steady states to different total copy numbers. Copy number variance initially increases sincetrajectories are all initialised at the same state, but plateaus because trajectories are constrained in theircopy number to remain near the attracting line of steady states (Figure S3A). Mean heteroplasmy remainsconstant through time under this model (Figure 2D, see (Birky et al. , 1983)). This is unsurprising since eachspecies possesses the same replication and degradation rate, so neither species is preferred.From stochastic simulations we observed that, for sufficiently short times, heteroplasmy variance increasesapproximately linearly through time for a range of parametrizations (Figure 2E-H), which is in agreementwith recent single-cell oocyte measurements in mice (Burgstaller et al. , 2018). Previous work has also shown alinear increase in heteroplasmy variance through time for purely genetic models of mtDNA dynamics (seeJohnston and Jones (2016)). We sought to understand the influence of mitochondrial network dynamics uponthe 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 ( x ( t = 0) = x ss), wetook the limit of limit of large mtDNA copy numbers, fast fission-fusion dynamics, and applied a second-ordertruncation of the Kramers-Moyal expansion (Gardiner, 1985) to the chemical master equation describingthe dynamics of the system (see Supporting Information). This yielded a stochastic differential equation for6 igure 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 steadystates under a deterministic model, as indicated by arrows. In stochastic simulation, mean copy number ( C ) is initiallyslightly perturbed from the deterministic treatment of the system, and then remains constant, while meanheteroplasmy ( D ) remains invariant with time (see Eq. (S61)). In ( E )-( H ), we show that Eq. (13) holds across manycellular circumstances: lines give analytic results, points are from stochastic simulation. Heteroplasmy variancebehaviour 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 γ → γ MR and β → β M where M and R denote the relative magnitude and ratio of the network rates, and γ , β denote the nominal parametrizations of thefusion and fission rates respectively (see Table S2). Figure S3D shows a sweep of M over the same logarithmic rangewhen R = 1. See Figure S4A-I and Table S3 for parameter sweeps numerically demonstrating the generality of theresult for different mtDNA control modes. et al. , 2016), we derived Eq. (S63), which may be approximated for sufficiently short times as V ( h ) ≈ f s ( x ) 2 µtn ( x ) h ( x )(1 − h ( x )) (cid:12)(cid:12)(cid:12)(cid:12) x = x ss . (13)Here, V ( h ) is the variance of heteroplasmy, µ is the mitophagy rate, n ( x ) is the total copy number and f s ( x ) is the fraction of unfused (singleton) mtDNAs, and is thus a measure of the fragmentation of themitochondrial network. x 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 fractionof unfused mitochondria, mitophagy rate, and inverse population size . We find that Eq. (13) closely matchesheteroplasmy variance dynamics from stochastic simulation, for sufficiently short times after initialisation, fora 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 thecellular population genetics of mtDNA variance. Its simple form allows several intuitive interpretations. Astime progresses, replication and degradation of both species occurs, allowing the ratio of species to fluctuate;hence we expect V ( h ) to increase with time according to random genetic drift (Figure 2E-H). The rate ofoccurrence of replication/degradation events is set by the mitophagy rate µ , since degradation events arebalanced by replication rates to maintain population size; hence, random genetic drift occurs more quickly ifthere is a larger turnover in the population (Figure 2E). We expect V ( h ) to increase more slowly in largepopulation 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 bythe 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 iszero. Furthermore, the rate of increase of heteroplasmy variance is maximal when a cell’s initial value ofheteroplasmy is 1/2. In Figure 2G, we show that Eq. (13) is able to recapitulate the rate of heteroplasmyvariance increase across different values of δ , which are hypothesized to correspond to different replicativesensing strengths of different mitochondrial mutations (Hoitzing, 2017). We also show in Figures S3B&C thatEq. (13) is robust to the choice of feedback control strength b in Eq. (10). n ( x ), f ( x ), and h ( x ) in Eq.(13)are not independent degrees of freedom in this model: they are functions of the state vector x , where x isdetermined by the parametrization and initial conditions of the model. Hence, the parameter sweeps in Figure2E-H and Figures S3B&C also implicitly vary over these functions of state by varying the steady state x ss.In Eq. (6), we have made the important assumption that only unfused mitochondria can be degradedvia mitophagy, as seen by Twig et al. (2008), hence the total propensity of mtDNA turnover is limited bythe number of mtDNAs which are actually susceptible to mitophagy. Strikingly, we find that the dynamicsof heteroplasmy variance are independent of the absolute rate of fusion and fission, only depending on thefraction 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 thatmitochondrial network dynamics are much faster than replication and degradation of mtDNA, by around afactor of β/µ ≈ (see Table S2), resulting in the existence of a separation of timescales between network andgenetic processes. In the derivation of Eq. (13), we have assumed that fission-fusion rates are infinite, whichsimplifies V ( h ) into a form which is independent of the magnitude of the fission-fusion rate. A parametersweep of the magnitude and ratio of the fission-fusion rates reveals that, if the fusion and fission rates aresufficiently small, Eq. (13) breaks down and V ( h ) gains dependence upon the magnitude of these rates (seeFigure S4A). This regime is, however, for network rates which are approximately 100 times smaller than thebiologically-motivated nominal parametrization shown in Figure 2A-D where the fission-fusion rate becomescomparable 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 ofeither 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 V ( h ) as heteroplasmy variance gradually becomes sub-linear with time, see Figure S5C&D. This is evidentthrough inspection of Eq. (S63), which shows that cellular trajectories which reach h = 0 or h = 1 cease todiffuse in heteroplasmy space, and so heteroplasmy variance cannot increase indefinitely. Consequently, thedepiction of heteroplasmy variance in Fig. 1B,D as being approximately normally distributed corresponds tothe regime in which our approximation holds, and is a valid subset of the behaviours displayed by heteroplasmydynamics under more sophisticated models (e.g. the Kimura distribution Kimura (1955); Wonnapinij et al. et al. (2008)). However, the linear regime for heteroplasmy variance has beenobserved 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 modelsof 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 ratefunction (i.e. λ = λ ( x ), µ = const), when the fusion and fission rates were high and the fixation probability( P ( h = 0) or P ( h = 1)) was negligible, Eq. (13) accurately described V ( h ) across all of the replicationrates investigated, see Figure S4A-F. The same mathematical argument to show Eq. (13) for the replicationrate in Eq. (10) may be applied to these alternative replication rates where a closed-form solution for thedeterministic steady state may be written down (see Supporting Information, Deriving an ODE description ofthe mitochondrial network system). Interestingly, when copy number is controlled through the degradationrate (i.e. λ =const, µ = µ ( x )), heteroplasmy variance loses its dependence upon network state entirely andthe f s term is lost from Eq. (13) (see Eq. (S72) and Figure S4G-I). A similar mathematical argument wasapplied to reveal how this dependence is lost (see Supporting Information, Proof of heteroplasmy relation forlinear feedback control).In order to provide an intuitive account for why control in the replication rate, versus control in thedegradation rate, determines whether or not heteroplasmy variance has network dependence, we investigateda time-rescaled form of the Moran process (see Supporting Information, A modified Moran process mayaccount for the alternative forms of heteroplasmy variance dynamics under different models of genetic mtDNAcontrol). The Moran process is structurally much simpler than the model presented above, to the point ofbeing unrealistic, in that the mitochondrial population size is constrained to be constant between consecutivetime steps. Despite this, the modified Moran process proved to be insightful. We find that, when copy numberis controlled through the replication rate, the absence of death in the fused subpopulation means the timescaleof the system (being the time to the next death event) is proportional to f s . In contrast, when copy number iscontrolled through the degradation rate, the presence of a constant birth rate in the entire population meansthe timescale of the system (being the time to the next birth event) is independent of f s (see Eq. (S84) andsurrounding discussion). Control strategies against mutant expansions
In this study, we have argued that the rate of increase of heteroplasmy variance, and therefore the rate ofaccumulation of pathologically mutated cells within a tissue, increases with mitophagy rate ( µ ), decreaseswith 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 theaccumulation of mutations. We use this new insight to propose three classes of strategy to control mutationaccumulation and hence address associated issues in aging and disease, and discuss these strategies throughthe 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, weinvoked an infinite sites Moran model (Kimura, 1969) (see Figure 3A). Single cells were modelled over timeas having a fixed mitochondrial copy number ( n ), and at each time step one mtDNA is randomly chosenfor 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 ∼ Binomial( L mtDNA , η ) (14)where Binomial( N, p ) is a binomial random variable with N trials and probability p of success. L mtDNA = 16569is the length of mtDNA in base pairs and η = 5 . × − 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 (Γ) of the Moran process is effectively rescaled by the fraction of unfused mitochondria,i.e. Γ = µnf s , which we apply here. 9 igure 3. Rate of de novo mutation accumulation is sensitive to the network state/mitophagy rateand 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 ) Meannumber of distinct mutations increases with the fraction of unfused mitochondria. This corresponds to a simplerescaling of time, so all but one of the parametrizations are shown in grey. ( C ) The mean number of mutations permtDNA also increases with the fraction of unfused mitochondria. Inset shows that the mean number of mutations permtDNA is independent of the number of mtDNAs per cell; values of n are the same as in ( D ). ( D ) Mean number ofmutations per cell increases according to the population size of mtDNAs. Standard error in the mean is too small tovisualise, so error bars are neglected, given 10 realizations. f s in Figures 3B&C therefore correspond to arescaling of time i.e. stretching of the time-axis. The absolute number of mutations predicted in Figure 3Bmay over-estimate the true number of mutations per cell (and of course depends on our choice of mutationrate), since a subset of mutations will experience either positive or negative selection. However, quantificationof the number of distinct mitochondrial mutants in single cells remains under-explored, as most mutationswill 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 inmitochondrial fusion (Mfn1 and Mfn2) in mouse skeletal muscle. Although knock-out studies presentdifficulties in extending their insights into the physiological case, the authors observed that fragmentationof the mitochondrial network induced severe depletion of mtDNA copy number (which we also observed inFigure S2B). Furthermore, the authors observed that the number of mutations per base pair increased uponfragmentation, which we also observed in the infinite sites model where fragmentation effectively results in afaster turnover of mtDNA (Figure 3C).Our models predict that promoting mitochondrial fusion has a two-fold effect: firstly, it slows the increaseof heteroplasmy variance (see Eq. (13) and Figure 2H); secondly, it reduces the rate of accumulation of distinctmutations (see Figure 3B&C). These two effects are both a consequence of mitochondrial fusion rescaling thetime to the next turnover event, and therefore the rate of random genetic drift. As a consequence, this simplemodel suggests that promoting fusion earlier in development (assuming mean heteroplasmy is low) could slowdown 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 leastfor some mutations under therapeutic conditions (Kandul et al. , 2016; Suen et al. , 2010), we predict thata 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 mitochondrialfusion and fission (Sukhorukov et al. , 2012; Zamponi et al. , 2018). In contrast, if selective mitophagy pathwaysare 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 µ have a comparable effect to changes in f s in terms of reducing the rateof heteroplasmy variance (see Eq. (13)) and the rate of de novo mutation (Figure 3B&C) since they bothserve to rescale time. Our theory therefore suggests that inhibition of basal mitophagy may be able to slowdown 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 stemcells (de Almeida et al. , 2017) have comparatively low levels of mitophagy, which is consistent with the ideathat these pluripotent cells attempt to minimise genetic drift by slowing down mtDNA turnover. A previousmodelling 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 healthywhen present at 100%, can induce fitness disadvantages (Acton et al. , 2007; Bagwan et al. , 2018; Sharpley et al. , 2012). In cases where heteroplasmy itself is disadvantageous, especially in later life where such mutationsmay have already accumulated, accelerating heteroplasmy variance increase to achieve fixation of a speciescould be advantageous. However, this will not avoid cell-to-cell variability, and the physiological consequencesfor 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 Γ = µn (i.e. a standard Moran process). We found that varying mtCN did not affect themean number of mutations per molecule of mtDNA (Figure 3C, inset). However, as the population sizebecomes larger, the total number of distinct mutations increases accordingly (Figure 3D). In contrast to ourpredictions, a recent study by Wachsmuth et al. (2016) found a negative correlation between mtCN and thenumber of distinct mutations in skeletal muscle. However, Wachsmuth et al. (2016) also found a correlation11etween the number of distinct mutations and age, in agreement with our model. Furthermore, the authorsused partial regression to find that age was more explanatory than mtCN in explaining the number of distinctmutations, 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 mtDNAmutations. Therefore, one would ideally account for these four variables for jointly, in order to fully constrainour model.A study of single neurons in the substantia nigra of healthy human individuals found that mtCN increasedwith age (D¨olle et al. , 2016). Furthermore, mice engineered to accumulate mtDNA deletions through faultymtDNA 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 theobserved 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 ina 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 coulddelay 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 anddegradation rates), it is often proposed that damaged mitochondria may experience a higher rate of degradation(Kim et al. , 2007; Narendra et al. , 2008). There are two principal ways in which selection may occur onmutant species. Firstly, mutant mitochondria may be excluded preferentially from the mitochondrial networkin a background of unbiased mitophagy. If this is the case, mutants would be unprotected from mitophagy forlonger periods of time than wild-types, and therefore be at greater hazard of degradation. We can alter thefusion rate ( γ ) in the mutant analogues of Eq. (1),(2) and Eqs. (7)–(9) by writing γ → γ/ (1 + (cid:15) f ) (15)for all fusion reactions involving 1 or more mutant mitochondria where (cid:15) f >
0. The second potential selectivemechanism we consider is selective mitophagy. In this case, the degradation rate of mutant mitochondria islarger than wild-types, i.e. we modify the mutant degradation reaction to M s µ (1+ (cid:15) m ) −−−−−→ ∅ (16)for (cid:15) m > (cid:15) f and (cid:15) m ) affects theextent of reduction in mean heteroplasmy. Figure 4A shows that, in the context of selective fusion ( (cid:15) f > )and non-selective mitophagy ( (cid:15) m = 0 ) the optimal strategy for clearance of mutants is to have an intermediatefusion/fission ratio . This was observed for all fusion selectivities investigated (see Figure S7) Intuitively, if themitochondrial network is completely fused then, due to mitophagy only acting upon smaller mitochondrialunits, mitophagy cannot occur – so mtDNA turnover ceases and heteroplasmy remains at its initial value. Incontrast, if the mitochondrial network completely fissions, there is no mitochondrial network to allow theexistence of a quality control mechanism: both mutants and wild-types possess the same probability per unittime of degradation, so mean heteroplasmy does not change. Since both extremes result in no clearance ofmutants, the optimal strategy must be to have an intermediate fusion/fission ratio.In contrast, in Figure 4B, in the context of non-selective fusion ( (cid:15) f = 0 ) and selective mitophagy ( (cid:15) m > ),the optimal strategy for clearance of mutants is to completely fission the mitochondrial network . Intuitively, ifmitophagy is selective, then the more mtDNAs which exist in fragmented organelles, the greater the numberof mtDNAs which are susceptible to selective mitophagy, the greater the total rate of selective mitophagy, thefaster the clearance of mutants. Discussion
In this work, we have sought to unify our understanding of three aspects of mitochondrial physiology – themitochondrial network state, mitophagy, and copy number – with genetic dynamics. The principal virtue ofour modelling approach is its simplified nature, which makes general, analytic, quantitative insights availablefor the first time. In using parsimonious models, we are able to make the first analytic link between the12 igure 4. Selective fusion implies intermediate fusion rates are optimal for mutant clearance whereasselective mitophagy implies complete fission is optimal . Numerical exploration of the shift in meanheteroplasmy for varying fusion/fission ratio, across different selectivity strengths. Stochastic simulations for meanheteroplasmy, evaluated at 1000 days, with an initial condition of h = 0 . n = 1000; the state was initialised on thesteady state line for the case of (cid:15) f = (cid:15) m = 0, for 10 iterations. ( A ) For selective fusion (see Eq. (15)), for each valueof fusion selectivity ( (cid:15) f ), the fusion rate ( γ ) was varied relative to the nominal parametrization (see Table S2). When (cid:15) f >
0, the largest reduction in mean heteroplasmy occurs at intermediate values of the fusion rate; a deterministictreatment reveals this to be true for all fusion selectivities investigated (see Figure S7). ( B ) For selective mitophagy(see Eq. (16)), when mitophagy selectivity (cid:15) 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. mitochondrial network state and heteroplasmy dynamics. This is in contrast to other computational studiesin the field, whose structural complexity make analytic progress difficult, and accounting for their predictedphenomena 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. Wefound, for a wide class of models of post-mitotic cells, that the rate of linear increase of heteroplasmy varianceis modulated in proportion to the fraction of unfused mitochondria (see Eq. (13)). The general notion thatmitochondrial fusion shields mtDNAs from turnover, and consequently serves to rescale time, emerges fromour analysis. This rescaling of time only holds when mitochondrial copy numbers are controlled througha state-dependent replication rate, and vanishes if copy numbers are controlled through a state-dependentmitophagy rate. We have presented the case of copy number control in the replication rate as being a moreintuitive model than control in the degradation rate. The former has the interpretation of biogenesis beingvaried to maintain a constant population size, with all mtDNAs possessing a characteristic lifetime. Thelatter has the interpretation of all mtDNA molecules being replicated with a constant probability per unittime, regardless of how large or small the population size is, and changes in mitophagy acting to regulatepopulation size. Such a control strategy seems wasteful in the case of stochastic fluctuations resulting in apopulation size which is too large, and potentially slow if fluctuations result in a population size which istoo small. Furthermore, control in the replication rate means that the mitochondrial network state may actas an additional axis for the cell to control heteroplasmy variance (Figure 2) and the rate of accumulationof de novo mutations (Figure 3B&C). Single-mtDNA tracking through confocal microscopy in conjunctionwith mild mtDNA depletion could shed light on whether the probability of degradation per unit time permtDNA varies when mtDNA copy number is perturbed, and therefore provide evidence for or against thesetwo possible control strategies.Our observations provide a substantial change in our understanding of mitochondrial genetics, as itsuggests that the mitochondrial network state, in addition to mitochondrial turnover and copy number, mustbe 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 ofheteroplasmy variance is independent of the absolute rate of fission-fusion events, since network dynamics13ccur approximately 10 times faster than mitochondrial turnover, inducing a separation of timescales. Theindependence of the absolute rate of network dynamics makes way for the possibility of gaining informationabout 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 rateof accumulation of de novo mutations, also due to the fraction of unfused mitochondria serving to rescaletime. We find that, in the context of mitochondrial quality control through selective fusion, an intermediatefusion/fission ratio is optimal due to the finite selectivity of fusion. This latter observation perhaps providesan indication for the reason why we observe mitochondrial networks in an intermediate fusion state underphysiological conditions (Sukhorukov et al. , 2012; Zamponi et al. , 2018).We have, broadly speaking, considered neutral models of mtDNA genetic dynamics. It is, however, typicallysuggested that increasing the rate of mitophagy promotes mtDNA quality control, and therefore shrinksthe distribution of heteroplasmies towards 0% mutant (see Eq. (15) and Eq. (16)). If mitophagy is able tochange mean heteroplasmy, then a neutral genetic model appears to be inappropriate, as mutants experiencea higher rate of degradation. Stimulation of the PINK1/Parkin pathway has been shown to select againstdeleterious mtDNA mutations in vitro (Suen et al. , 2010) and in vivo (Kandul et al. , 2016), as has repressionof the mTOR pathway via treatment with rapamycin (Dai et al. , 2013; Kandul et al. , 2016). However, thenecessity of performing a genetic/pharmacological intervention to clear mutations via this pathway suggeststhat 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 heteroplasmyevolves through time under physiological conditions. Indeed, it has been recently shown that mitophagy isbasal (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 importantfor understanding the role of mitochondrial mutations in healthy aging (Kauppila et al. , 2017; Khrapko andVijg, 2009). A typical rate of increase of heteroplasmy variance predicted by Eq.(13) given our nominalparametrization (Table S2) is V (cid:48) ( h ) /t = V ( h ) / ( E ( h )(1 − E ( h )) t ) = 2 µf s /n ≈ . × − day − ( f s = 0 . n = 1000). This value accounts for the accumulation of heteroplasmy variance which is attributable to turnoverof the mitochondrial population in a post-mitotic cell. However, in the most general case, cell division is alsoable to induce substantial heteroplasmy variance. For example, V (cid:48) ( h ) /t has been measured in model organismgermlines to be approximately 9 × − day − in Drosophila (Johnston and Jones, 2016; Solignac et al. ,1987), 9 × − day − in NZB/BALB mice (Johnston and Jones, 2016; Wai et al. , 2008; Wonnapinij et al. ,2008), and 2 × − day − in single LE and HB mouse oocytes (Burgstaller et al. , 2018). We see that theserates of increase in heteroplasmy variance are approximately an order of magnitude larger than predictionsfrom our model of purely quiescent turnover, given our nominal parametrisation. Whilst larger mitophagyrates may also potentially induce larger values for V (cid:48) ( h ) /t (see Poovathingal et al. (2012), and Figure S5C,corrsponding to V (cid:48) ( h ) /t ≈ . × − day − ) 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 (Aryaman et al. , 2019; Kauppila et al. , 2017).Our findings reveal some apparent differences with previous studies which link mitochondrial geneticswith network dynamics (see Table S4). Firstly, Tam et al. (2013, 2015) found that slower fission-fusiondynamics resulted in larger increases in heteroplasmy variance with time, in contrast to Eq. (13) which onlydepends 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 multiplemtDNAs. Faster fission-fusion dynamics tended to form heteroplasmic mitochondria whereas slower dynamicsformed homoplasmic mitochondria. It is intuitive that mitophagy of a homoplasmic mitochondrion induces alarger shift in heteroplasmy than mitophagy of a single mtDNA, hence slower network dynamics form morehomoplasmic mitochondria. However, this apparent difference with our findings can naturally be resolved if weconsider 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 (Burgstaller et al. , 2014a; Cagalinec et al. , 2013). If the fission-fusionrates are sufficiently large to ensure heteroplasmic mitochondria, then further increasing the fission-fusion rateis unlikely to have an impact on heteroplasmy dynamics. Hence, this finding is potentially compatible withour study, although future experimental studies investigating intra-mitochondrial heteroplasmy would helpconstrain these models. Tam et al. (2015) also found that fast fission-fusion rates could induce an increasein mean heteroplasmy, in contrast to Figure 2D which shows that mean heteroplasmy is constant with time14fter a small initial transient due to stochastic effects. We may speculate that the key difference between ourtreatment and that of Tam et al. (2013, 2015) is the inclusion of cellular subcompartments which inducesspatial effects which we do not consider here. The uncertainty in accounting for the phenomena observed insuch complex models highlights the virtues of a simplified approach which may yield interpretable laws andprinciples through analytic treatment.The study of Mouli et al. (2009) suggested that, in the context of selective fusion, higher fusion ratesare optimal. This initially seems to contrast with our finding which states that intermediate fusion ratesare optimal for the clearance of mutants (Figure 4A). However, the high fusion rates in that study do notcorrespond 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 optimalclearing are an order of magnitude lower than the highest fusion rate they consider. In the case of completefusion, mitophagy cannot occur in the model of Mouli et al. (2009), so there is no mechanism to removedysfunctional 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 selectivemitophagy 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 fusionrates 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 maynot be functional; ii) the sigmoidal relationship between number of functional units per mitochondrion andmitochondrial ‘activity’ (the metric by which optimality is measured). Points (i) and (ii) imply that smallnumbers of dysfunctional mitochondrial units have very little impact on mitochondrial activity, so fusion mayboost total mitochondrial activity in the context of small amounts of mutation. So whilst Figure 4B remainsplausible in light of the study of Mouli et al. (2009) if reduction of mean heteroplasmy is the objective of thecell, 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 andselective mitophagy. Future experimental studies quantifying the importance of selective mitophagy underphysiological conditions would be beneficial for understanding heteroplasmy variance dynamics. The ubiquityof heteroplasmy (Morris et al. , 2017; Payne et al. , 2012; Ye et al. , 2014) suggests that a neutral drift approachto mitochondrial genetics may be justified, which contrasts with the studies of Tam et al. (2013, 2015) andMouli 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 studyby 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, arerequired to account for the rate of increase in heteroplasmy variance. Joint longitudinal measurements of f s , µ and n , with heteroplasmy quantification, would allow verification of Eq. (13) and aid in determining the extentto which neutral genetic models are explanatory. This could be achieved, for instance, using the mito-QCmouse (McWilliams et al. , 2016) which allows visualisation of mitophagy and mitochondrial architecture in vivo . Measurement of f s , µ and n , followed by e.g. destructive single-cell whole-genome sequencing ofmtDNA would allow validation of how µ , n and f s influence V ( h ) and the rate of de novo mutation (seeFigure 3). One difficulty is sequencing errors induced through e.g. PCR, which hampers our ability accuratelymeasure 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 mitochondrionpossessing 3.9 single-nucleotide variants on average. Error correction strategies during sequencing may pavethe way towards high-accuracy mtDNA sequencing (Salk et al. , 2018; Woods et al. , 2018), and allow us tobetter constrain models of heteroplasmy dynamics. Acknowledgments
We would like to thank Hanne Hoitzing, Thomas McGrath, Abhishek Deshpande, and Ferdinando Insalatafor useful discussions. Simulations were performed using the Imperial College High Performance ComputingService. J.A. acknowledges grant support from the BBSRC (BB/J014575/1), and the MRC MitochondrialBiology Unit (MC UP 1501/2). N.S.J. acknowledges grant support from the BHF (RE/13/2/30182) andEPSRC (EP/N014529/1). I.G.J. acknowledges funding from ERC StG 805046 (EvoConBiO) and a TuringFellowship from the Alan Turing Institute. 15 eferences
Acton, B., I. Lai, X. Shang, A. Jurisicova, and R. Casper, 2007 Neutral mitochondrial heteroplasmy altersphysiological function in mice. Biol. Reprod. : 569–576.Aryaman, J., I. G. Johnston, and N. S. Jones, 2017 Mitochondrial DNA density homeostasis accounts for athreshold effect in a cybrid model of a human mitochondrial disease. Biochem. J. : 4019–4034.Aryaman, J., I. G. Johnston, and N. S. Jones, 2019 Mitochondrial heterogeneity. Front. Genet. : 718.Assaf, M. and B. Meerson, 2010 Extinction of metastable stochastic populations. Phys. Rev. E : 021116.Bagwan, N., E. Bonzon-Kulichenko, E. Calvo, A. V. Lechuga-Vieco, S. Michalakopoulos, et al. , 2018Comprehensive quantification of the modified proteome reveals oxidative heart damage in mitochondrialheteroplasmy. Cell Rep. : 3685–3697.Birky, C. W., T. Maruyama, and P. Fuerst, 1983 An approach to population and evolutionary genetic theoryfor genes in mitochondria and chloroplasts, and some results. Genetics : 513–527.Boudoures, A. L., J. Saben, A. Drury, S. Scheaffer, Z. Modi, et al. , 2017 Obesity-exposed oocytes accumulateand transmit damaged mitochondria due to an inability to activate mitophagy. Dev. Biol. : 126–138.Burgstaller, J. P., I. G. Johnston, N. S. Jones, J. Albrechtov´a, T. Kolbe, et al. , 2014a MtDNA segregation inheteroplasmic tissues is common in vivo and modulated by haplotype differences and developmental stage.Cell Rep. : 2031–2041.Burgstaller, J. P., I. G. Johnston, and J. Poulton, 2014b Mitochondrial DNA disease and developmentalimplications for reproductive strategies. Mol. Hum. Reprod. : 11–22.Burgstaller, J. P., T. Kolbe, V. Havlicek, S. Hembach, J. Poulton, et al. , 2018 Large-scale genetic analysisreveals mammalian mtDNA heteroplasmy dynamics and variance increase through lifetimes and generations.Nat. Commun. : 2488.Cagalinec, M., D. Safiulina, M. Liiv, J. Liiv, V. Choubey, et al. , 2013 Principles of the mitochondrial fusionand fission cycle in neurons. J. Cell Sci. : 2187–2197.Capps, G. J., D. C. Samuels, and P. F. Chinnery, 2003 A model of the nuclear control of mitochondrial DNAreplication. J. Theor. Biol. : 565–583.Chen, H., M. Vermulst, Y. E. Wang, A. Chomyn, T. A. Prolla, et al. , 2010 Mitochondrial fusion is requiredfor mtDNA stability in skeletal muscle and tolerance of mtDNA mutations. Cell : 280–289.Chinnery, P. F. and D. C. Samuels, 1999 Relaxed replication of mtDNA: a model with implications for theexpression of disease. Am. J. Hum. Genet. : 1158.Constable, G. W., T. Rogers, A. J. McKane, and C. E. Tarnita, 2016 Demographic noise can reverse thedirection of deterministic selection. Proc. Natl. Acad. Sci. USA p. 201603693.Cree, L., S. Patel, A. Pyle, S. Lynn, D. Turnbull, et al. , 2008 Age-related decline in mitochondrial DNA copynumber in isolated human pancreatic islets. Diabetologia : 1440–1443.Dai, Y., K. Zheng, J. Clark, R. H. Swerdlow, S. M. Pulst, et al. , 2013 Rapamycin drives selection against apathogenic heteroplasmic mitochondrial DNA mutation. Hum. Mol. Genet. : 637–647.de Almeida, M. J., L. L. Luchsinger, D. J. Corrigan, L. J. Williams, and H.-W. Snoeck, 2017 Dye-independentmethods reveal elevated mitochondrial mass in hematopoietic stem cells. Cell Stem Cell : 725–729.D¨olle, C., I. Flønes, G. S. Nido, H. Miletic, N. Osuagwu, et al. , 2016 Defective mitochondrial DNA homeostasisin the substantia nigra in Parkinson disease. Nat. Commun. : 13548.Durham, S. E., D. C. Samuels, L. M. Cree, and P. F. Chinnery, 2007 Normal levels of wild-type mitochondrialDNA maintain cytochrome c oxidase activity for two pathogenic mitochondrial DNA mutations but not form.3243A → G. Am. J. Hum. Genet. : 189–195. 16lson, J., D. Samuels, D. Turnbull, and P. Chinnery, 2001 Random intracellular drift explains the clonalexpansion of mitochondrial dna mutations with age. Am. J. Hum. Genet. : 802–806.Figge, M. T., A. S. Reichert, M. Meyer-Hermann, and H. D. Osiewacz, 2012 Deceleration of fusion–fissioncycles improves mitochondrial quality control during aging. PLoS Comp. Biol. : e1002576.Gardiner, C., 1985 Stochastic Methods: A Handbook for the Natural and Social Sciences . Springer Series inSynergetics, Springer Berlin Heidelberg.Gillespie, D. T., 1976 A general method for numerically simulating the stochastic time evolution of coupledchemical reactions. J. Comp. Phys. : 403–434.Gillespie, D. T., 2007 Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem. : 35–55.Grima, R., 2010 An effective rate equation approach to reaction kinetics in small volumes: Theory andapplication to biochemical reactions in nonequilibrium steady-state conditions. J. Chem. Phys. : 07B604.Grima, R., P. Thomas, and A. V. Straube, 2011 How accurate are the nonlinear chemical Fokker-Planck andchemical Langevin equations? J. Chem. Phys. : 084103.Hoitzing, H., 2017 Controlling mitochondrial dynamics: population genetics and networks . Ph.D. thesis,Imperial College London.Hoitzing, H., P. A. Gammage, M. Minczuk, I. G. Johnston, and N. S. Jones, 2017 Energetic costs of cellularand therapeutic control of stochastic mtDNA populations. bioRxiv p. 145292.Hoitzing, H., I. G. Johnston, and N. S. Jones, 2015 What is the function of mitochondrial networks? Atheoretical assessment of hypotheses and proposal for future research. BioEssays : 687–700.Jacobs, K., 2010 Stochastic processes for physicists: understanding noisy systems . Cambridge UniversityPress.Johnston, I. G., 2018 Tension and resolution: Dynamic, evolving populations of organelle genomes withinplant cells. Mol. Plant .Johnston, I. G., J. P. Burgstaller, V. Havlicek, T. Kolbe, T. R¨ulicke, et al. , 2015 Stochastic modelling, bayesianinference, and new in vivo measurements elucidate the debated mtDNA bottleneck mechanism. eLife :e07464.Johnston, I. G. and N. S. Jones, 2016 Evolution of cell-to-cell variability in stochastic, controlled, heteroplasmicmtDNA populations. Am. J. Hum. Genet. : 1150–1162.Kandul, N. P., T. Zhang, B. A. Hay, and M. Guo, 2016 Selective removal of deletion-bearing mitochondrialDNA in heteroplasmic drosophila. Nat. Commun. : 13100.Kauppila, T. E., J. H. Kauppila, and N.-G. Larsson, 2017 Mammalian mitochondria and aging: an update.Cell Metab. : 57–71.Khrapko, K. and J. Vijg, 2009 Mitochondrial DNA mutations and aging: devils in the details? Trends Genet. : 91–98.Kim, I., S. Rodriguez-Enriquez, and J. J. Lemasters, 2007 Selective degradation of mitochondria by mitophagy.Arch. Biochem. Biophys. : 245–253.Kimura, M., 1955 Solution of a process of random genetic drift with a continuous model. Proc. Natl. Acad.Sci. USA : 144.Kimura, M., 1969 The number of heterozygous nucleotide sites maintained in a finite population due to steadyflux of mutations. Genetics : 893–903.Kowald, A. and T. Kirkwood, 2018 Resolving the enigma of the clonal expansion of mtdna deletions. Genes :126. 17owald, A. and T. B. Kirkwood, 2013 Mitochondrial mutations and aging: random drift is insufficient toexplain the accumulation of mitochondrial deletion mutants in short-lived animals. Aging Cell : 728–731.Kowald, A. and T. B. Kirkwood, 2014 Transcription could be the key to the selection advantage of mitochondrialdeletion mutants in aging. Proc. Natl. Acad. Sci. USA : 2972–2977.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, C., C. A. Wurm, H. Sp˚ahr, M. Falkenberg, N.-G. Larsson, et al. , 2011 Super-resolution microscopyreveals that mammalian mitochondrial nucleoids have a uniform size and frequently contain a single copyof mtDNA. Proc. Natl. Acad. Sci. USA : 13534–13539.Lewis, S. C., L. F. Uchiyama, and J. Nunnari, 2016 ER-mitochondria contacts couple mtDNA synthesis withmitochondrial division in human cells. Science : aaf5549.Li, M., R. Schr¨oder, S. Ni, B. Madea, and M. Stoneking, 2015 Extensive tissue-related and allele-relatedmtDNA heteroplasmy suggests positive selection for somatic mutations. Proc. Natl. Acad. Sci. USA :2491–2496.Lieber, T., S. P. Jeedigunta, J. M. Palozzi, R. Lehmann, and T. R. Hurd, 2019 Mitochondrial fragmentationdrives selective removal of deleterious mtDNA in the germline. Nature p. 1.L´opez-Ot´ın, C., M. A. Blasco, L. Partridge, M. Serrano, and G. Kroemer, 2013 The hallmarks of aging. Cell : 1194–1217.McWilliams, T. G., A. R. Prescott, G. F. Allen, J. Tamjar, M. J. Munson, et al. , 2016 mito-QC illuminatesmitophagy and mitochondrial architecture in vivo. J. Cell Biol. pp. jcb–201603039.McWilliams, T. G., A. R. Prescott, L. Montava-Garriga, G. Ball, F. Singh, et al. , 2018 Basal mitophagyoccurs independently of PINK1 in mouse tissues of high metabolic demand. Cell Metab. : 439–449.Medeiros, T. C., R. L. Thomas, R. Ghillebert, and M. Graef, 2018 Autophagy balances mtDNA synthesis anddegradation by DNA polymerase POLG during starvation. J. Cell Biol. pp. jcb–201801168.Mengel-From, J., M. Thinggaard, C. Dalg˚ard, K. O. Kyvik, K. Christensen, et al. , 2014 Mitochondrial DNAcopy number in peripheral blood cells declines with age and is associated with general health among elderly.Hum. Genet. : 1149–1159.Mishra, P. and D. C. Chan, 2016 Metabolic regulation of mitochondrial dynamics. J. Cell Biol. : 379–387.Morris, J., Y.-J. Na, H. Zhu, J.-H. Lee, H. Giang, et al. , 2017 Pervasive within-mitochondrion single-nucleotidevariant heteroplasmy as revealed by single-mitochondrion sequencing. Cell Rep. : 2706–2713.Mouli, P. K., G. Twig, and O. S. Shirihai, 2009 Frequency and selectivity of mitochondrial fusion are key toits quality maintenance function. Biophys. J. : 3509–3518.Narendra, D., A. Tanaka, D.-F. Suen, and R. J. Youle, 2008 Parkin is recruited selectively to impairedmitochondria and promotes their autophagy. J. Cell Biol. : 795–803.Pan, J., L. Wang, C. Lu, Y. Zhu, Z. Min, et al. , 2019 Matching mitochondrial DNA haplotypes for circumventingtissue-specific segregation bias. iScience : 371–379.Parsons, T. L. and T. Rogers, 2017 Dimension reduction for stochastic dynamical systems forced onto amanifold by large drift: a constructive approach with examples from theoretical biology. J. Phys. A :415601.Patel, P. K., O. Shirihai, and K. C. Huang, 2013 Optimal dynamics for quality control in spatially distributedmitochondrial networks. PLoS Comp. Biol. : e1003108.Payne, B. A., I. J. Wilson, P. Yu-Wai-Man, J. Coxhead, D. Deehan, et al. , 2012 Universal heteroplasmy ofhuman mitochondrial DNA. Hum. Mol. Genet. : 384–390.18erier, C., A. Bender, E. Garc´ıa-Arum´ı, M. J. Melia, J. Bove, et al. , 2013 Accumulation of mitochondrialDNA deletions within dopaminergic neurons triggers neuroprotective mechanisms. Brain : 2369–2378.Picard, M., J. Zhang, S. Hancock, O. Derbeneva, R. Golhar, et al. , 2014 Progressive increase in mtDNA3243A > G heteroplasmy causes abrupt transcriptional reprogramming. Proc. Natl. Acad. Sci. USA :E4033–E4042.Poovathingal, S. K., J. Gruber, B. Halliwell, and R. Gunawan, 2009 Stochastic drift in mitochondrial DNApoint mutations: a novel perspective ex silico. PLoS Comput. Biol. : e1000572.Poovathingal, S. K., J. Gruber, L. Lakshmanan, B. Halliwell, and R. Gunawan, 2012 Is mitochondrial DNAturnover slower than commonly assumed? Biogerontology : 557–564.Rossignol, R., B. Faustin, C. Rocher, M. Malgat, J. Mazat, et al. , 2003 Mitochondrial threshold effects.Biochem. J. : 751–762.Salk, J. J., M. W. Schmitt, and L. A. Loeb, 2018 Enhancing the accuracy of next-generation sequencing fordetecting rare and subclonal mutations. Nature Rev. Genet. .Samuels, D. C., C. Li, B. Li, Z. Song, E. Torstenson, et al. , 2013 Recurrent tissue-specific mtDNA mutationsare common in humans. PLoS Genet. : e1003929.Samuels, D. C., P. Wonnapinij, L. Cree, and P. F. Chinnery, 2010 Reassessing evidence for a postnatalmitochondrial genetic bottleneck. Nat. Genet. : 471–472.Schon, E. A., S. DiMauro, and M. Hirano, 2012 Human mitochondrial DNA: roles of inherited and somaticmutations. Nat. Rev. Genet. : 878–890.Sharpley, M. S., C. Marciniak, K. Eckel-Mahan, M. McManus, M. Crimi, et al. , 2012 Heteroplasmy of mousemtDNA is genetically unstable and results in altered behavior and cognition. Cell : 333–343.Solignac, M., J. G´enermont, M. Monnerot, and J.-C. Mounolou, 1987 Drosophila mitochondrial genetics:evolution of heteroplasmy through germ line cell divisions. Genetics : 687–696.Stewart, J. B. and P. F. Chinnery, 2015 The dynamics of mitochondrial DNA heteroplasmy: implications forhuman health and disease. Nature Rev. Genet. : 530–542.Stewart, J. B. and N.-G. Larsson, 2014 Keeping mtDNA in shape between generations. PLoS Genet. :e1004670.Suen, D.-F., D. P. Narendra, A. Tanaka, G. Manfredi, and R. J. Youle, 2010 Parkin overexpression selectsagainst a deleterious mtDNA mutation in heteroplasmic cybrid cells. Proc. Natl. Acad. Sci. USA :11835–11840.Sukhorukov, V. M., D. Dikov, A. S. Reichert, and M. Meyer-Hermann, 2012 Emergence of the mitochondrialreticulum from fission and fusion dynamics. PLoS Comput. Biol. : e1002745.Szabadkai, G., A. M. Simoni, K. Bianchi, D. De Stefani, S. Leo, et al. , 2006 Mitochondrial dynamics andCa2+ signaling. Biochim. Biophys. Acta. : 442–449.Tam, Z. Y., J. Gruber, B. Halliwell, and R. Gunawan, 2013 Mathematical modeling of the role of mitochondrialfusion and fission in mitochondrial DNA maintenance. PloS ONE : e76230.Tam, Z. Y., J. Gruber, B. Halliwell, and R. Gunawan, 2015 Context-dependent role of mitochondrialfusion-fission in clonal expansion of mtDNA mutations. PLoS Comp. Biol. : e1004183.Taylor, R. W., M. J. Barron, G. M. Borthwick, A. Gospel, P. F. Chinnery, et al. , 2003 Mitochondrial DNAmutations in human colonic crypt stem cells. J. Clin. Invest. : 1351–1360.Trifunovic, A., A. Wredenberg, M. Falkenberg, J. N. Spelbrink, A. T. Rovio, et al. , 2004 Premature ageing inmice expressing defective mitochondrial DNA polymerase. Nature : 417.Twig, G., A. Elorza, A. J. Molina, H. Mohamed, J. D. Wikstrom, et al. , 2008 Fission and selective fusiongovern mitochondrial segregation and elimination by autophagy. EMBO J. : 433–446.19an Kampen, N. G., 1992 Stochastic processes in physics and chemistry , volume 1. Elsevier.Wachsmuth, M., A. Huebner, M. Li, B. Madea, and M. Stoneking, 2016 Age-related and heteroplasmy-relatedvariation in human mtDNA copy number. PLoS Genet. : e1005939.Wai, T., D. Teoli, and E. A. Shoubridge, 2008 The mitochondrial DNA genetic bottleneck results fromreplication of a subpopulation of genomes. Nat. Genet. : 1484.Wallace, D. C. and D. Chalkia, 2013 Mitochondrial DNA genetics and the heteroplasmy conundrum inevolution and disease. Cold Spring Harb. Perspect. Biol. : a021220.Wilkinson, D. J., 2011 Stochastic modelling for systems biology . CRC press.Wonnapinij, P., P. Chinnery, and D. Samuels, 2008 The distribution of mitochondrial DNA heteroplasmy dueto random genetic drift. Am. J. Hum. Genet. : 582.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 : 265.Ye, K., J. Lu, F. Ma, A. Keinan, and Z. Gu, 2014 Extensive pathogenicity of mitochondrial heteroplasmy inhealthy human individuals. Proc. Natl. Acad. Sci. USA : 10654–10659.Zamponi, N., E. Zamponi, S. A. Cannas, O. V. Billoni, P. R. Helguera, et al. , 2018 Mitochondrial networkcomplexity emerges from fission/fusion dynamics. Sci. Rep. : 363.Zheng, W., K. Khrapko, H. A. Coller, W. G. Thilly, and W. C. Copeland, 2006 Origins of human mitochondrialpoint mutations as DNA polymerase γ -mediated errors. Mutat. Res. : 11–20.20 upporting InformationConstant rates yield unstable copy numbers for a model describingmtDNA genetic and network dynamics We explored a simpler network system than the one presented in the Main Text, but found that it producedinstability in mtDNA copy numbers, which we regard as biologically undesirable. Consider the following setof Poisson processes for singleton ( s ) and fused ( f ) species s + s γ −→ f + f (S1) s + f γ −→ f + f (S2) f β −→ s (S3) s αρ −→ s + s (S4) f αρ −→ f + f (S5) s ηρ −→ ∅ (S6) f ρ −→ ∅ (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 αρ . Eq. (S6) and (S7) are simple death processeswith rates ηρ and ρ respectively. The parameter ρ is shared amongst all of the birth and death reactions inEqs. (S4)–(S7). ρ represents the intuitive assumption that, in order for a stable population size to exist, birthshould balance death. However, for the network to have any effect at all, singletons should be at an increasedrisk of mitophagy relative to fused species. We represent the increased risk of singleton mitophagy with theparameter η . Since additional death is introduced into the system when η >
1, we include the parameter α > s d t = − γs − γf s + βf + αρs − ηρs (S8)d f d t = γs + γf s − βf + αρf − ρf (S9)where we have enforced the stochastic reaction rate to be equivalent to the deterministic reaction rate, andhence the s term is proportional to γ rather than 2 γ (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 saddlenode, 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 examplesuggests that a system of this form with constant reaction rates is unstable, and therefore biologically unlikelyto exist under reasonable circumstances. We hence consider analogous biochemical reaction networks with areplication 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 differentialequations
The following section outlines the steps in converting a set of chemical reactions into a set of ordinary differentialequations (ODEs). In particular, we pay special attention to the fact that the rate of a chemical reaction witha 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 anumber of standard texts (Gillespie, 1976, 2007; Van Kampen, 1992; 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 j X + · · · + s Nj X N ˆ k j −→ r j X + · · · + r Nj X N (S10)21here s ij and r ij are stoichiometric coefficients. We define ˆ k j as the microscopic rate for this reaction. Thedimensionality of this parameter will vary depending upon the stoichiometric coefficients s ij . ˆ k j may beloosely 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 thesystem and time, moving forwards through time. Defining the state of the system as x = ( x , . . . , x N ) T ,where x i is the copy number of the i th species, allows us to write the CME as (Grima, 2010) ∂P ( x , t | x , t ) ∂t = Ω R (cid:88) j =1 (cid:32) N (cid:89) i =1 E − S ij i − (cid:33) ˆ f j ( x , Ω) P ( x , t | x , t ) (S11)where Ω is the volume of the compartment in which the reactions occur (also known as the system size), S ij = r ij − s ij is the stoichiometry matrix, and E − S ij i is referred to as the step operator and is definedthrough the relation E − S ij i ( g ( x )) = g ( x , . . . , x i − S ij , . . . , x N ), for any function of state g ( x ). ˆ f j ( x , Ω) is themicroscopic rate function of reaction j , which in general depends on both the state and the system size. Afactor of Ω is explicitly included in this definition of the chemical master equation so that our treatment iscompatible with Van Kampen’s system size expansion (Van Kampen, 1992). As a consequence of this, theprobability that, given the current state x , the j th reaction occurs in the time interval [ t, t + d t ) somewhere inΩ (Gillespie, 2007) is ˆ a j ( x , Ω)d t := Ω ˆ f j ( x , Ω)d t. (S12)ˆ a j ( x , Ω) is termed the propensity function (or “hazard”) and is of particular relevance in the stochasticsimulation algorithm (Gillespie, 1976), since ˆ a j ( x , Ω) / (cid:80) j ˆ a j ( x , Ω) determines the probability that the j th reaction occurs next.For the microscopic rate function, we may writeˆ f j ( x , Ω) = ˆ k j N (cid:89) i =1 Ω − s ij (cid:18) x i s ij (cid:19) . (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 ofthe CME) which is defined as (Grima, 2010; Van Kampen, 1992)d φ i d t = R (cid:88) j =1 S ij ˜ f j ( φ ) (S14)where φ = ( φ , . . . , φ N ) T is the vector of macroscopic concentrations (of dimensions molecules per unitvolume) and ˜ f j ( φ ) is the macroscopic rate function satisfying˜ f j ( φ ) = ˜ k j N (cid:89) i =1 φ s ij i (S15)where ˜ k j is the macroscopic rate for the j th reaction. We distinguish between ˆ k j and ˜ k j , respectively the rateconstants for the discrete and continuous pictures, although this distinction is sometimes not emphasizedin the literature (Grima, 2010; Grima et al. , 2011; Van Kampen, 1992). The physical meaning of ˜ k j isnot immediately obvious: we argue that this parameter only gains physical meaning through the followingprocedure.As stated by Wilkinson (2011), if we intend for the microscopic description in Eq. (S11) to correspond tothe macroscopic description in Eq. (S14), the rate of consumption/production of particles for every reactionmust be the same in the deterministic limit of the stochastic system (the conditions for which we definebelow). Therefore, we apply the following constraint in the limit of large copy numberslim x i →∞ ˆ f j ( x , Ω) = ˜ f j ( φ ) ∀ i, j. (S16)22n applying this constraint on all species i and all reactions j , we may derive a general relationship betweenˆ k j and ˜ k j lim x i →∞ ˆ k j N (cid:89) i =1 Ω − s ij x i ! s ij ! ( x i − s ij )! = ˜ k j N (cid:89) i =1 φ s ij i . (S17)We can make two approximations to generate a more convenient relationship between the microscopic andmacroscopic rates. Firstly, we assume that x i ≈ Ω φ i . (S18)This is a small noise approximation, since it is often assumed that x i = Ω φ i + Ω / ξ i , where ξ i is a noise term(Van Kampen, 1992). If ξ i is small then x i ≈ Ω φ i is a valid approximation. Secondly, we assume that x i ( x i − . . . ( x i − s ij + 1) ≈ x s ij i . (S19)This is a large copy number approximation: in the case of e.g. a bimolecular reaction (2 X i → ∗ ) with s ij = 2,the approximation is of the form x i ( x i − ≈ x i or x i ≈ x i −
1. By applying Eq. (S19) to the factor of x i ! inEq. (S17), the factor of ( x i − s ij )! cancels from the left-hand side. Simplifying using Eq. (S18), φ s ij i cancelsfrom both sides and we arrive at the important relationship˜ k j ≈ ˆ k j N (cid:89) i =1 s ij ! . (S20)With Eq. (S14), Eq. (S15) and Eq. (S20) one may therefore write down a set of ODEs for an arbitrary chemicalreaction network, with constant reaction rates, in terms of the microscopic rates ˆ k j . This equation highlightsthat for reactions with s ij ≥
2, ˜ k j (cid:54) = ˆ k j , as is the case for bimolecular reactions of the form 2 X i → ∗ (seeEq. (1) and Eq. (S1)).Importantly, if the microscopic rate function is a function of state then ˆ k = ˆ k ( x ) and ˜ k = ˜ k ( φ ) ≈ ˜ k ( x / Ω).In this case, Eq. (S20) still applies since the above argument assumed nothing about the particular forms of ˆ k and ˜ k . However, additional factors of Ω − are induced by applying Eq. (S18), which may carry through tothe individual parameters of ˜ k ( φ ). A demonstration of this is given in the following section. Deriving an ODE description of the mitochondrial network system
In this section we show how to derive an ODE description of the network system described in Eq. (1)-Eq. (9)in the Main Text. In accordance with the notation in the previous section, we will redefine all of the rates inEq. (1)-Eq. (10) with a hat notation (ˆ a , for a general rate parameter a ), to reflect that these are stochasticrates. Deterministic rates will be denoted with a tilde (˜ a ). Our aim will be to write a set of ODEs in terms ofthe stochastic rates, ˆ 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 rewriteEq. (1) to allow the reaction to proceed with some arbitrary rate ˆ ρ : X S + X S ˆ ρ −→ X F + X F , (S21)where X denotes either a wild-type ( W ) or mutant ( M ). We will subsequently fix ˆ ρ to the rate of all otherfusion reactions ˆ γ . We do this because Eq. (1) is a bimolecular reaction involving one species: a fundamentallydifferent reaction to bimolecular reactions involving two species, as we will now see.Since ˆ ρ, ˆ γ = const, we may use Eq. (S20), resulting in the deterministic rates˜ ρ = ˆ ρ γ = ˆ γ (S23)for Eq. (S21) and Eq. (2) respectively. If we then enforce the microscopic rates to be equal for both of thesefusion reactions, i.e. ˆ ρ = ˆ γ , then ˜ ρ = ˆ γ/
2. All other fusion reactions have ˜ γ = ˆ γ by application of Eq. (S20).Application of Eq. (S20) to the fission reaction in Eq. (3) shows that ˆ β = ˜ β .23or Eq. (4), we have chosen a ˆ λ which is not a constant, but a function of the copy numbers of the chemicalspecies (ˆ λ = ˆ λ ( x ) where x = ( w s , w f , m s , m f ), see Eq. (10)). As pointed out in the previous section, caremust be taken in writing down the deterministic analogue of ˆ λ . Applying Eq. (S20), we haveˆ µ + ˆ b (ˆ κ − ( w s + w f + ˆ δm s + ˆ δm f )) = ˜ µ + ˜ b (˜ κ − ( φ w s + φ w f + ˜ δφ m s + ˜ δφ m f )) . (S24)Applying Eq. (S18) and equating individual terms, we arrive atˆ µ = ˜ µ (S25)ˆ bw s = ˜ b w s Ω = ⇒ ˆ b = ˜ b/ Ω (S26)ˆ b ˆ k = ˜ b ˜ k = ⇒ ˆ k = ˜ k Ω (S27)ˆ b ˆ δm s = ˜ b ˜ δφ m s = ⇒ ˆ δ = ˜ δ. (S28)In this study, we let Ω = 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 ODEsd φ w s d t = − · ˆ γ φ w s − ˆ γφ w s φ w f + ˆ βφ w f − (ˆ µ + ˆ b (ˆ κ − ( φ w s + φ w f +ˆ δφ m s + ˆ δφ m f ))) φ w s − ˆ µφ w s − ˆ γφ m f φ w s − ˆ γφ w s φ m s (S29)d φ m s d t = − · ˆ γ φ m s − ˆ γφ m s φ m f + ˆ βφ m f − (ˆ µ + ˆ b (ˆ κ − ( φ w s + φ w f +ˆ δφ m s + ˆ δφ m f ))) φ m s − ˆ µφ m s − ˆ γφ w f φ m s − ˆ γφ w s φ m s (S30)d φ w f d t = 2 · ˆ γ φ w s + ˆ γφ w s φ w f − ˆ βφ w f + (ˆ µ + ˆ b (ˆ κ − ( φ w s + φ w f +ˆ δφ m s + ˆ δφ m f )))(2 φ w s + φ w f )+ˆ γφ m f φ w s + ˆ γφ w s φ m s (S31)d φ m f d t = 2 · ˆ γ φ m s + ˆ γφ m s φ m f − ˆ βφ m f + (ˆ µ + ˆ b (ˆ κ − ( φ w s + φ w f +ˆ δφ m s + ˆ δφ m f )))(2 φ m s + φ m f )+ˆ γφ w f φ m s + ˆ γφ w s φ m s . (S32)The steady state solution of this system of ODEs may be calculated, but its form is complex. For notationalsimplicity, we will drop the hat notation. Defining x = ( b ( β + 2 β ( γκ + 3 µ ) + γ κ + µ + 2 γµ ( κ + 2( δ − φ m s )) + 2 bγµ ( − β + µ + γ ( κ − δφ m s + 2 φ m s )) + γ µ ) / (S33)the non-trivial, physically-realizable, component of the steady state may be parametrized in terms of φ m s andwritten as φ w s = − ( βb κ + b γκ + b κµ + 2 b δµφ m s + β b + βbγκ + 5 βbµ + 3 bγκµ − βγµ + 2 bµ − bγδµφ m s − bγµφ m s − x ( bκ + β + 2 µ ) + 2 γµ + 2 γ µφ m s ) / (2 µ ( b − γ ) ) (S34)24 w f = ( b ( β κ + β ( γκ ( κ + ( δ − φ m s ) + µ (3 κ + δφ m s )) + φ m s ( γ ( δ − κ − δµ + γµ ((2 δ − κ + 2( δ − δφ m s ))) + b ( β + β (6 µ + γ ( κ + ( δ − φ m s )) + β (5 µ + γ δκφ m s + γ κ ( − φ m s ) + 6 γδµφ m s − γµφ m s − κx )+ φ m s (2 γ µ ( δκ + δ ( − φ m s ) + φ m s ) + γ ((6 δ − µ − ( δ − κx ) − δµx )) − b ( β ( γ µ ( − ( κ + (4 − δ ) φ m s )) + γ (2 µ + ( δ − φ m s x ) + 3 µx ) + γµφ m s ( γ ( δ − µ + γ ( κ − δ − φ m s ) + ( δ − x )+ β (2 γµ + x )) + γµ ( β − γφ m s )( γµ + x )) / (2 bµ ( b − γ ) ( β + γ ( δ − φ m s )) (S35) φ m f = ( φ m s ( − bβ + bγκ + bµ − bγδφ m s + 2 bγφ m s + γµ + x )) / (2 b ( β + γ ( δ − φ m s )) . (S36)Since the steady state is parametrized by φ 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 replicationrate in Eq. (10) using the Kramers-Moyal expansion under conditions of large copy number and fast networkchurn (to be defined below); the approach used here is similar to Constable et al. (2016). Consonant withthe self-contained objectives of STAR methods, we draw together elements from the literature to provide acoherent 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 initiallyproceed in this way. Following the treatment by Gardiner (1985), we begin by re-writing the chemical masterequation Eq. (S11) (CME) as ∂P ( x , t ) ∂t = (cid:90) ∞−∞ d x (cid:48) [ T ( x | x (cid:48) ) P ( x (cid:48) , t ) − T ( x (cid:48) | x ) P ( x , t )] (S37)where we have set Ω = 1. T ( x | x (cid:48) ) is the transition rate from state x (cid:48) → x , and the dependence upon theinitial condition has been suppressed for notational convenience. We now proceed by expanding the CME.The multivariate Kramers-Moyal expansion may be written as ∂P ( x , t ) ∂t ≈ (cid:90) ∞−∞ (cid:18) −∇ ( T ( x (cid:48) | x ) P ( x )) T · ( x (cid:48) − x ) + 12 ( x (cid:48) − x ) T · H · ( x (cid:48) − x ) (cid:19) d x (cid:48) (S38)where H ( x ) is the Hessian matrix of T ( x (cid:48) | x ) P ( x ) H := ∂ ∂x . . . ∂ ∂x ∂x N ... ... ∂ ∂x N ∂x . . . ∂ ∂x N T ( x (cid:48) | x ) P ( x ) (S39)(see (Gardiner, 1985) for a proof of this in the univariate case).A transition to each possible neighbouring state x (cid:48) corresponds to some reaction j which moves the statefrom x → x (cid:48) . Since we know the influence of each reaction on state x through the constant stoichiometrymatrix S ij , and that the propensity of a reaction does not depend upon x (cid:48) itself (see Eq. (S13)), we maytransition from a notation involving x and x (cid:48) into a notation involving x and j . We may therefore define T j ( x ) := T ( x (cid:48) | x ) ≡ ˆ f j ( x ) (see Eq. (S13)), and let H ( x ) → H j ( x ).25e now make a large copy number assumption in order to simplify T j ( x ). To take a large copy numberlimit, we assume that x i ! ≈ x s ij i ( x i − s ij )! resulting in T j ( x ) ≈ ˆ k j N (cid:89) i =1 x s ij i s ij ! . (S40)This approximation is exact when s ij = 0 ,
1, but inexact when s ij ≥
2. For example, if we consider the second-order bimolecular reaction in Eq. (1), Eq. (S40) is equivalent to assuming w s ≈ w s ( w s − / ( s ij ! ) = 1 / T j ( x ) as a combinatorial factor from stochastic considerations. Fokker-Planck equation for chemical reaction networks
We now wish to re-write Eq. (S38) as aFokker-Planck equation. Since the integral in Eq. (S38) is over x (cid:48) , and every x (cid:48) corresponds to a reaction j ,we may interpret the integral in Eq. (S38) as a sum over all reactions, i.e. (cid:82) d x (cid:48) → (cid:80) Rj =1 . Hence, for the j th reaction, [( x (cid:48) − x )] i = S ij . With these observations, we may write the first integral of Eq. (S38) as (cid:90) ∞−∞ −∇ ( T ( x (cid:48) | x ) P ( x )) T · ( x (cid:48) − x ) d x (cid:48) = (cid:90) ∞−∞ −∇ ( T j ( x ) P ( x , t )) T · ( x (cid:48) − x ) d x (cid:48) = − R (cid:88) j =1 N (cid:88) i =1 ∂∂x i ( T j ( x ) P ( x , t )) S ij = − N (cid:88) i =1 ∂∂x i A i P ( x , t ) (S41)where A := S · T . (S42) A is a vector of length N , [ S ] ij := r ij − s ij is the N × R stoichiometry matrix Eq. (S10), and T is the vectorof 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 H j in Eq. (S39) as H jlm = ∂ ∂x l ∂x m T j ( x ) P ( x , t ) (S43)where j = 1 , . . . , R and l, m = 1 , . . . , N . Thus, we may write (cid:90) ∞−∞
12 ( x (cid:48) − x ) T · H j · ( x (cid:48) − x ) d x (cid:48) = 12 R (cid:88) j =1 N (cid:88) l =1 N (cid:88) m =1 S lj H jlm S mj = 12 R (cid:88) j =1 N (cid:88) l =1 N (cid:88) m =1 S lj ∂ ∂x l ∂x m T j P ( x , t ) S mj = 12 N (cid:88) l =1 N (cid:88) m =1 ∂ ∂x l ∂x m R (cid:88) j =1 S lj T j S mj P ( x , t )= 12 N (cid:88) i,m =1 ∂ ∂x i ∂x m B im P ( x , t ) (S44)where B := S · Diag( T ) · S T . (S45) B is an N × N matrix, and Diag( Y ) is a diagonal matrix whose main diagonal is the vector Y . We maytherefore re-write Eq. (S38) as a Fokker-Planck equation for the state vector x of the form ∂P ( x , t ) ∂t ≈ − N (cid:88) i =1 ∂∂x i [ A i ( x ) P ( x , t )] + 12 N (cid:88) i,m =1 ∂ ∂x i ∂x m [ B im ( x ) P ( x )] . (S46)26 okker-Planck equation for an arbitrary function of state We now wish to make a change of variablesin Eq. (S46) to write down a Fokker-Planck equation for an arbitrary scalar function of state x (which we willlater set to be heteroplasmy). To do this, we wish to make use of Itˆo’s formula, which allows a change ofvariables for an SDE. In general, the Fokker-Planck equation in Eq. (S46) is equivalent (Jacobs, 2010) to thefollowing Itˆo stochastic differential equation (SDE)d x = A d t + G d W (S47)where GG T ≡ B (where G is an N × R matrix) and d W is a vector of independent Wiener increments oflength R , and a Wiener increment d W satisfies (cid:90) t d W := W ( t ) , P ( W, t ) ≡ √ πt e − W / (2 t ) . (S48)Itˆo’s formula states that, for an arbitrary function h ( x , t ) where x satisfies Eq. (S47), we may write thefollowing SDE d h ( x , t ) = (cid:26) ∂h∂t + ( ∇ h ) T A + 12 Tr (cid:2) G T H h ( x ) G (cid:3)(cid:27) d t + ( ∇ h ) T G d W , (S49)where H h ( x ) is the Hessian matrix of h ( x , t ) (see Eq. (S39), where T ( x (cid:48) | x ) P ( x ) should be replaced with h ( x , t )). Given the form of B in Eq. (S45) we let G = S · Diag( √ T ) , (S50)which satisfies GG T ≡ B .For convenience, we may also perform the transformation purely at the level of Fokker-Planck equations.Let h ( x , t ) satisfy the general Fokker-Planck equation ∂P ( h, t ) ∂t = − ∂∂h [ ˜ A ( h, t ) P ( h, t )] + 12 ∂ ∂h [ ˜ B ( h, t ) P ( h, t )] (S51)for scalar functions ˜ A ( h, t ) and ˜ B ( h, t ). Using the cyclic property of the trace in Eq. (S49), we may identify˜ A = ∂h∂t + ( ∇ h ) T A + 12 Tr [ B H h ( x )] (S52)where Tr is the trace operator. Also, from Eq. (S49),˜ B = [( ∇ h ) T G ][( ∇ h ) T G ] T = ( ∇ h ) T B ( ∇ h ) . (S53)Hence, using Eq. (S51), Eq. (S52) and Eq. (S53), we may write down a Fokker-Planck equation for an arbitraryfunction of state in terms of A and B . An SDE for heteroplasmy forced onto the steady state line in the high-churn limit
It has beendemonstrated that SDE descriptions of stochastic systems which possess a globally-attracting line of steadystates 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 tracesout 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 traceout the position on the steady state line. We seek to use similar reasoning to verify Eq. (13). In what follows,we will assume that x ( t = 0) = x ss, where x ss is the state which is the solution of A = (which is equivalentto finding the steady state solution of the deterministic rate equation in Eq. (S14) due to our assumption oflarge copy numbers and Ω = 1), so that we may neglect any deterministic transient dynamics.Inspection of the steady state of the ODE description of our system reveals that the set of steady statesolutions forms a line (see Eqs. (S34)–(S36)). Inspection of the steady state solution reveals that the steadystate depends on the fusion ( γ ) and fission ( β ) rates. Mitochondrial network dynamics occur on a muchfaster timescale than the replication and degradation of mtDNA: the former occurring on the timescale ofminutes (Twig et al. , 2008) whereas the latter is hours or days (Johnston and Jones, 2016). We seek to use27his separation of timescales to arrive at a simple form for V ( h ). We redefine the fusion and fission rates suchthat γ → M γβ → M β (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 ( x , t ) = h ( x ) := ( 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 thediffusion term ˜ B for heteroplasmy using Eq. (S53). If we constrain the state x to be forced onto the steadystate line x ss (as per (Constable et al. , 2016; Parsons and Rogers, 2017)) in the high-churn limit, then upondefining θ := (cid:112) b ( β + γκ ) − bγµ ( β − γκ + 2 γ ( δ − m s ) + γ µ (S56)we have lim M →∞ (cid:18) ˜ B (cid:12)(cid:12)(cid:12) x = x ss (cid:19) = (16 b γ µm s ( β + γ ( δ − m s ) ( b ( β + 2 β γ ( δ − m s + βγ ( κ + (1 − δ ) κm s + ( δ − δ − m s ) + γ κm s (( δ − m s − κ )) + b ( − β (2 γµ + θ ) + βγ ( κ (2 γµ + θ ) + m s ( γ (5 − δ ) µ − δ − θ )) + γ m s (( δ − m s (3 γµ + θ ) − κ (2 γµ + θ ))) + γµ ( γµ + θ )( β − γm s ))) / ( β ( b ( γ ( κ − δm s +2 m s ) − β ) + γµ + θ ) ) . (S57)Eq. (S57) is difficult to understand. In order to perform further simplification, we make an ansatz for theform of ˜ B ( ˜ B An ) and seek to determine whether our ansatz is equivalent to the derived form of ˜ B under theconstraints defined on the left-hand side of Eq. (S57). Our ansatz takes the form˜ B An := lim M →∞ (cid:32) µh (1 − h ) n ( x ) · f s ( x ) (cid:12)(cid:12)(cid:12)(cid:12) x = x ss (cid:33) (S58)where f s ( x ) := ( w s + m s ) / ( w s + w f + m s + m f ) and n ( x ) := w s + w f + m s + m f . Notice that this ansatzis more general than Eq. (S57), since it has no explicit dependence upon the parameters of the control lawassumed in Eq. (10), and only explicitly depends upon functions of state x .Upon substituting the steady state x ss into the ansatz in Eq. (S58) and taking the high-churn limit, wefind that ˜ B An = − ( m s ( β + γ ( δ − m s )( b ( β + γκ ) − γµ + θ )( b ( β + γκ ) + γµ − θ )( b ( m s (2 βδ − β + γκ ) − βκ ) + m s ( θ − γµ ))) / (2 bβ ( κ − δm s + m s ) ( b ( γ ( κ − δm s + 2 m s ) − β ) + γµ + θ )) . (S59)After some algebra (see GitHub repository for Mathematica notebook), it can be shown that Eq. (S57) andEq. (S59) are equivalent, i.e. ˜ B An ≡ lim M →∞ (cid:18) ˜ B (cid:12)(cid:12)(cid:12) x = x ss (cid:19) . (S60)28s such, we may use ˜ B and ˜ B An interchangeably in the limit of high network churn. Furthermore, it can beshown after some algebra that the drift of heteroplasmy when forced onto the steady state line is 0, i.e.˜ A (cid:12)(cid:12)(cid:12) x = x ss ≡ . (S61)A similar result is shown in (Constable et al. , 2016) ((Equation S59) therein). Substituting h ( x , t ) = h , ˜ A and˜ B into the Fokker-Planck equation for an arbitrary function of state Eq. (S51), we have ∂P ( h, t ) ∂t = 12 ∂ ∂h (cid:34) (cid:18) µh (1 − h ) n ( x ) · f s ( x ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) x = x ss( h ) P ( h, t ) (cid:35) (S62)which is equivalent to the following SDE for heteroplasmyd h = (cid:115) µh (1 − h ) f s ( x ) n ( x ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = x ss( h ) d W (S63)in the limit of large network churn, large copy numbers, and a second-order truncation of the Kramers-Moyalexpansion. Although the state has been forced onto the steady state, stochastic fluctuations mean thattrajectories may move along the line of steady states, so the diffusion coefficient is not constant in general.We may calculate the new value of x ss( h ) for every displacement due to Wiener noise in h , and substituteinto f s ( x ) and n ( x ) to determine the diffusion coefficient at the next time step.However, for sufficiently short times, and large copy numbers (i.e. low diffusivity of h ), we may assumethat the diffusion coefficient in Eq. (S63) may be approximated as constant. Since the general solution of theSDE d y = √ B d W (S64)for B = const is y ∼ N ( y | y , Bt ) (S65)where N ( y | y , σ ) is a Gaussian distribution on y with mean y and variance σ , and y = y ( t = 0). Since wehave assumed that the state is initialised at x ( t = 0) = x ss, there are no deterministic transient dynamics, sowe may write V ( h ) ≈ µtn ( x ) h ( x )(1 − h ( x )) f s ( x ) (cid:12)(cid:12)(cid:12)(cid:12) x = x ss , (S66)where V returns variance of a random variable. In this equation, we take x = x ss = const, since we haveassumed a low-diffusion limit. We observe that this equation is of precisely the same form as (Equation 12) ofJohnston and Jones (2016), except with an additional proportionality factor of f s induced by the inclusion ofa mitochondrial network. Heteroplasmy variance relations for alternative model structuresand modes of genetic mtDNA control
Here we explore the implications of alternative model structures upon Eq. (S63). Firstly, we may considerreplacing Eq. (4) with X S λ −→ X S + X S . (S67)This corresponds to the case where replication coincides with fission, see (Lewis et al. , 2016). Repeating thecalculation in the previous section also results in Eq. (S63), so the result is robust to the particular choice ofmtDNA replication reaction (see GitHub repository for Mathematica notebook).Secondly, we may explore the impact of allowing non-zero mtDNA degradation of fused species. Thiscould correspond to autophagy-independent degradation of mtDNA, for example via the exonuclease activityof POLG (Medeiros et al. , 2018). To encode this, we may add the following additional reaction X F ξµ −→ ∅ (S68)29here 0 ≤ ξ ≤
1. We were not able to make analogous analytical progress in this instance. However, numericalinvestigation (Figure S3E) revealed that the following ansatz was able to predict heteroplasmy variancedynamics V ( h ) ≈ µtn ( x ) h ( x )(1 − h ( x ))( f s ( x ) + ξ (1 − f s ( x ))) (cid:12)(cid:12)(cid:12)(cid:12) x = x ss . (S69)In other words, allowing degradation of fused species results in a linear correction to our heteroplasmy varianceformula in Eq. (13). If fused species are susceptible to degradation at the same rate as unfused species ( ξ = 1),then V ( h ) loses f s dependence entirely and the mitochondrial network has no influence over heteroplasmydynamics.We also explored various different forms of λ ( x ) and µ ( x ), which we label A-G after (Johnston and Jones,2016), and X-Z for several newly-considered functional forms, see Table 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 discussionsurrounding Eq. (S1). The argument presented in the previous section requires the steady state solution ofthe system to be solvable, since we require the explicit form of x 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 theabove can be applied (see the GitHub repository for Mathematica notebooks). Controls B, C, E, F all satisfyEq. (S63); this can be shown numerically for controls A and X. However, controls G, Y and Z satisfyd h = (cid:115) µh (1 − h ) n ( x ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = x ss( h ) 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 intuitionfor why the diffusion rate of heteroplasmy variance depends on the network state when the population iscontrolled through replication, and does not depend on network state when the population is controlledthrough 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 ± β = 33 .
12 day − .The dimensions of β are day − and not mitochondrion − day − . This is because if the propensity (seeEq. (S12), where Ω = 1) of e.g. Eq. (3) is ˆ a fis ,w = βw f then the mean time to the next event is 1 / ( βw f );therefore the dimension of β is per unit time and copy numbers are pure numbers, i.e. dimensionless. Similarreasoning constrains the dimension of the fusion rate, see below.Evaluation of the fusion rate is more involved, since fusion involves two different chemical species comingtogether to react whereas fission may be considered as spontaneous. Furthermore, there are 7 differentfusion reactions whereas there are only 2 fission reactions. For simplicity, assume that all species have asteady-state copy number of x i = 250 (resulting in a total copy number of 1000, heteroplasmy of 0.5 and50% of mitochondria existing in the fused state). Neglecting subtleties relating to bimolecular reactionsinvolving one species (see Eq. (S20)), each fusion reaction proceeds at rate ˆ a fus ,j ≈ γx i . Since there are7 fusion reactions (Eq. (1), Eq. (2), Eq. (7)-Eq. (9)), the total fusion propensity is ˆ a fus ≈ γx i . Similarly,the total fission propensity is ˆ a fis = β ( w f + m f ) = 2 βx i . Since we expect macroscopic proportions of bothfused and fissioned species in many physiological settings, we may equate the fusion and fission propensities,ˆ a fus = ˆ a fis , and rearrange for the fusion rate γ to yield γ = 2 βx i / (7 x i ) ≈ . × − day − . The orders ofmagnitude difference between β and γ stems from the observation that fusion propensity depends on thesquare of copy number whereas the fission propensity depends on copy number linearly.Given the network parameters, we then explored appropriate parametrizations for the genetic parameters:the mitophagy rate ( µ ) and the parameters of the linear feedback control ( κ , b and δ , see Eq. (10)). mtDNAhalf-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 andheteroplasmy, we use an mtDNA half-life of 30 days (Tam et al. , 2015).30he parameter δ in the replication feedback control (see Eq. (10)) may be interpreted as the “strength ofsensing of mutant mtDNA” in the feedback control (Hoitzing et al. , 2017). Assuming that fluctuations in copynumber of mutants and wild-type molecules are sensed identically (as may be the case for e.g. non-codingmtDNA mutations) we may reasonably assume a model of δ = 1 as the simplest case of a neutral mutation(although δ (cid:54) = 1 still defines a neutral model, since both mutant and wild-type alleles experience the samereplication and degradation rates per molecule, see Eqs.(4)–(6)).We are finally left with setting the parameters κ and b in the linear feedback control Eq. (10). In theabsence of a network state and mutants, κ is precisely equal to the steady state copy number, since thedegradation rate equals the replication rate when κ = w . However, the presence of a network means that asubpopulation of mtDNAs (namely the fused species) are immune to death, resulting in κ no longer beingequivalent to the steady state copy number. The parameter b may be interpreted as the feedback controlstrength, which determines the extent to which the replication rate changes given a unit change in copynumber.Given a particular value of b , we may search for a κ which gives a total steady state copy number ( n )which is closest to some target value (e.g. 1000 as a typical total mtDNA copy number per cell in humanfibroblasts (Kukat et al. , 2011)). We swept a range of different values of b and found that, for values of b smaller than a critical value ( b ∗ ), a κ could not be found whose deterministic steady state was sufficientlyclose to n = 1000. This result is intuitive because in the limit of b → λ = const. From the analysis abovewe have shown that constant genetic rates ( µ, λ ) result in unstable copy numbers, and therefore a sufficientlysmall value of b is not expected to yield a stable non-trivial steady state solution. We chose b ≈ b ∗ , and thecorresponding κ , such that the steady state copy number is controlled as weakly as possible given the modelstructure. Rate renormalization
In Eqs. (1)–(10) we have neglected reactions such as X F + X F → X F + X F (S71)because they do not change the number of molecules in our state vector x = ( w s , w f , m s , m f ). One may askwhether neglecting such reactions means that it is necessary to renormalize the fission-fusion rates which wereestimated in the preceding section. In estimating the nominal parametrization above, we began by using aliterature value for the mitochondrial fission rate, and then matched the fusion rate such that the summedhazard of a fusion event approximately balanced the fission rate. This matching procedure is reasonable, sincewe observe a mixture of fused and fissioned mitochondria under physiological conditions: choosing a fusionrate which is vastly different results in either a hyperfused or fragmented network. We must therefore onlyjustify the fission rate. Eq. (3) assumes that a fission reaction always results in a singleton, and a singletonis by definition a molecule which is susceptible to mitophagy (see Eq. (6)). Therefore, if fission reactionsalways result in mitochondria containing single mtDNAs which are susceptible to mitophagy, then we expectour model to match well to true physiological rates. If, on the other hand, fission reactions between largecomponents of the network which are too large to be degraded are common, then renormalization of β by thefraction of fission events which result in a sufficiently small mitochondrion would be necessary, which would inturn renormalize γ through our matching procedure. We are not aware of experimental measurements of thefraction of fission events which result in mitochondria which contain a particular number of mtDNAs. Suchan experiment, combined with the distribution of mitochondrial sizes which are susceptible to mitophagy,would allow us to validate our approach. Despite this, the robustness of our results over approximately 4orders of magnitude for the fission-fusion rate (Figure S4A-I) provides some indication that our results arelikely to hold in physiological regions of parameter space. A modified Moran process may account for the alternative forms ofheteroplasmy variance dynamics under different models of geneticmtDNA control
We sought to gain insight into why control of population size through the replication rate ( λ = λ ( x ), µ =const) results in heteroplasmy variance depending on the fraction of unfused mitochondria (see Eq. (13)),31hereas control of population size through the degradation rate ( µ = µ ( x ), λ = const) results in heteroplasmyvariance becoming independent of network state, where V ( h ) ≈ λtn ( x ) h ( x )(1 − h ( x )) (cid:12)(cid:12)(cid:12)(cid:12) x = x ss . (S72)We will proceed by considering an analogous Moran process to the set of reactions presented in Eqs. (1) (9).First, consider a haploid biallelic Moran process consisting of wild-types and mutants, in a population offixed size n . At each step in discrete time, a member of the population is chosen for birth, and another fordeath. Let m t denote the copies of mutants at time t . Then, P ( m t +1 = j | m t = i ) = i ( n − i ) /n if j = i ± i /n + ( n − i ) /n if j = i E ( m t +1 | m t ) = m t . (S74)Defining h t := m t /n then from Eq. (S73) V ( m t +1 | m t ) = 2 h t (1 − h t ) (S75)and therefore V (cid:16) m t +1 n (cid:12)(cid:12)(cid:12) m t (cid:17) = V ( h t +1 | m t ) = 2 n h t (1 − h t ) . (S76)Suppose that, instead of the process occurring with discrete time, instead the process occurs with continuoustime, where each event is a simultaneous birth and death, and is modelled as a Poisson process. Suppose thatevents occur at a rate µ per capita. The waiting time between successive events ( τ ) is an exponential randomvariable with rate µN . Hence the expected waiting time between successive events is E ( τ ) = 1 µn . (S77)If we take the ratio of Eq. (S76) and Eq. (S77), we have V ( h t + τ | m t ) E ( τ ) = 2 µh t (1 − h t ) n . (S78)Heuristically, one could interpret Eq. (S78) as a ratio of differentials as follows. If we were to suppose that n were large enough such that E ( τ ) is very small, and h t is approximately constant ( h ) after a small number ofevents, then ∆ V ( h )∆ τ ≈ V ( h t + τ | m t ) E ( τ ) = ⇒ d V ( h )d t ≈ µh (1 − h ) n = ⇒ V ( h, t ) ≈ µh (1 − h ) n t (S79)where we have replaced the inter-event time τ with physical time t . This result is analogous to Eq. (S72) andEq. (12) of Johnston and Jones (2016), and agrees with simulation (Figure S6A).Now consider the modified Moran process in Figure S6B, which we refer to as a “protected” Moran process.Let 0 < f s ≤ nf s h t and nf s (1 − h t )mutants and wild-types are randomly chosen to be susceptible to death, respectively, where n is large. In thiscontinuous-time model, the inter-event time is τ ∼ Exponential(Γ) where Γ will be defined below. Then anindividual from the susceptible population is chosen for death, and any individual is allowed to be born. Thebirth and death events occur simultaneously in time.Again, using t as an integer counter of events, we have P ( m t +1 = j | m t = i ) = nf s h t nf s n (1 − h t ) N if j = i − nf s (1 − h t ) nf s nh t N if j = i + 1 nf s h t nf s nh t N + nf s (1 − h t ) nf s n (1 − h t ) n if j = i h t (1 − h t ) if j = i ± h t + (1 − h t ) if j = i τ ∼ Exponential(Γ). Firstly, if the death rate per capita is constant ( µ ), then the rate at which a death eventoccurs in the system (Γ death ) is proportional to the number of individuals which are susceptible to death:Γ death = µnf s . If we assume that the overall birth rate is matched to the overall death rate so that populationsize is maintained, as is the case when λ = λ ( x ) in the network system, then the overall birth rate (Γ birth )must also be Γ birth = µnf s . Hence, Γ = Γ birth + Γ death = 2 µnf s (S82)where µ is a proportionality constant. Since, for a Moran event to occur, both a birth and a death event mustoccur, time effectively runs twice as fast in a Moran model relative to a comparable chemical reaction networkmodel. We therefore rescale time by taking µ → µ/
2, and thusΓ = µnf s . (S83)As a result, E ( τ ) = 1 / ( µnf s ) and therefore, using Eq. (S76) and the reasoning in Eq. (S79),∆ V ( h )∆ τ ≈ V ( h t + τ | m t ) E ( τ ) = ⇒ d V ( h )d t ≈ µf s h (1 − h ) n = ⇒ V ( h ) ≈ µf s h (1 − h ) n t. (S84)This is analogous to when λ = λ ( x ) and µ = const in the network system. Hence, when λ = λ ( x ) and µ = const, the absence of death in the fused subpopulation means the timescale of the system (being the time to thenext death event) is proportional to f s . This argument is only a heuristic, since the Moran process is definedsuch that birth and death events occur simultaneously and therefore do not possess separate propensities(Γ birth and Γ death ).The second case we consider is when each individual has a constant rate of birth, hence Γ birth ∝ n . Thenthe death rate is chosen such that Γ birth = Γ death . In this case Γ = λn , where λ is a proportionality constant.The same argument from Eq. (S77) to Eq. (S79) may be applied, with an appropriate rescaling of time, and wearrive at Eq. (S72). This is analogous to when µ = µ ( x ) and λ = const in the network system. Hence, when µ = µ ( x ) and λ = const, the presence of a constant birth rate in the entire population means the timescale ofthe system (being the time to the next birth event) is independent of f s .33 able S1. Key predictions from our mathematical models . The following results hold for our neutral genetic model of a post-mitotic cell, with a simple model ofmitochondrial network dynamics:1. The rate of increase of heteroplasmy variance is proportional to the fraction of unfused mitochondria,but independent of the absolute magnitude of fission-fusion rates, due to a rescaling of time by themitochondrial network (Eq.(13), Eq.(S63), Figure 2E-H).2. The rate of accumulation of de novo mutations increases as the fraction of unfused mitochondriaincreases, due to a rescaling of time by the mitochondrial network (Figure 3B-D)3. When fusion is selective, intermediate fusion-fission ratios are optimal for reducing mean heteroplasmy(Figure 4A)4. When mitophagy is selective, complete fission is optimal for reducing mean heteroplasmy (Figure 4B)
Figure S1. Phase portrait for an ODE representation of a network system with constant rates . Thesystem displays two steady states: a trivial steady state at s + f = 0, and a non-trivial steady state. At a steady state,both time derivatives in Eq. (S8) and Eq. (S9) vanish. Trajectories (blue) show the evolution of the system up to t = 1000. The direction and magnitude of the derivative at points in space are shown by red arrows. Trajectories canbe seen either decaying to s = f = 0 or tending to infinity. γ = 0 . β = 33 . ρ = 0 . η = 1 . α = 1 . able S2. Nominal parametrization for network system . Dimensions of parameters are derived usingindividual terms in Eqs. (S29)–(S32); copy numbers of particular species are pure numbers and thereforedimensionless. Since we have chosen a system size of Ω = 1 throughout, we set the dimension of volume to 1. SeeChoice of nominal parametrization for further details of parameter justification. Parameter Description Value Dimensions Remarks β Fission rate 33.12 day − For cortical neurons, see(Cagalinec et al. , 2013) γ Fusion rate 3 . × − day − Approximately balancesfission µ Mitophagy rate 0.023 day − From (Tam et al. , 2015) δ Mutant feedback sensi-tivity 1 dimensionless Potentially appropriatefor e.g. non-codingmtDNA mutations b Feedback controlstrength 1 . × − day − Chosen as the weakestcontrol strength whichhas a non-trivial steadystate and total copynumber of 1000 κ Steady state copy num-ber parameter 11.7 dimensionless See remark for b Figure S2. Deterministic treatment of network system . ( A ) Deterministic dynamics of total copy numberunder linear feedback control, which is controlled to a particular steady state value (see Figure 2A&B). ( B ) Defining aknock-down (KD) factor ( k − = 0 . , . . . . , . β → β/k (red) and the fusion rate to γ → γ/k (blue), causing a linear increase and decrease in total copy number respectively under a deterministictreatment (see Figure 2A&B). igure S3. Stochastic treatment of network system . ( A ) Copy number variance for stochastic simulationsinitially increases, since all stochastic simulations begin with the same initial condition, but then plateaus since thesteady state line is globally attracting (see Figure 2C&D). ( B ) Error in Eq. (13) in a sweep over the feedback controlstrength, b . Dotted line denotes a 5% error according to Eq. (12). ( C ) V ( h ) profile for the parametrization with thelargest error in ( B ). ( D ) Sweeps of the network rate magnitude (see Figure 2H). Heteroplasmy variance isapproximately independent of absolute network rates over a broad range of network magnitudes. ( E ) Allowing fusedspecies to be degraded with relative rate ξ (Eq. (S68)), stochastic simulations for heteroplasmy variance (markers) andEq. (S69) (lines) are shown. Fused species degradation induces a linear correction to the heteroplasmy varianceformula. igure S4. Parameter sweeps of network fission-fusion rates for replication-based anddegradation-based control modes show robustness of two heteroplasmy variance formulaerespectively . ( A-I ) Error in Eq. (13) (ansatz) and Eq. (S72) (ansatz indep network) for sweeps of the fusion andfission rates for the corresponding feedback control functions in Table S3. Equations are accurate to at least 5% (blueregions) across large regions of parameter space, for many control laws. Fusion and fission rates are redefined as γ → γ MR and β → β M where M and R denote the magnitude and ratio of the network rates, and γ , β denotethe nominal parametrizations of the fusion and fission rates respectively. Summary statistics for 10 realizations, withinitial condition h = 0 . t = 500 days. Errors in V ( h ) (see Eq. (12)) smaller than 5% are truncatedand are shown as blue. Parametrizations where a deterministic steady state could not be found for an initial conditionof h = 0 . h = 0. Where insets arenot present, the probability of fixation is negligible. ( A-F ) When λ = λ ( x ) and P ( h = 0) is low, Eq. (13) performs wellin the high-churn limit. ( G-I ) When µ = µ ( x ) and P ( h = 0) is low, Eq. (S72) performs well in the high-churn limit. able S3. Nominal parametrizations for the alternative feedback control functions explored . Nominalparametrizations for the feedback controls in Fig. S4. In all cases, the nominal fission and fusion rates were β = 33 . γ = 0 .
038 respectively.
Interpretation Replication rate Degradation rate Note
Linear feedback (see Fig-ure 2, S2, S3, S4A) µ + b ( κ − w T − δm T ); seeTable S2 µ ; see Table S2 Control E in (Johnstonand Jones, 2016) andGitHub repositoryRelaxed replication (seeFigure S4B) αµ ( w opt − w T − δm T ) / ( w T + m T ); α = 1, w opt = 1000, δ = 1 µ ; µ = 0 .
023 Control A in (Johnstonand Jones, 2016) andGitHub repositoryDifferential control fortarget population (seeFigure S4C) α ( w opt − w T ) µ ; µ = 0 .
023 Control B in (Johnstonand Jones, 2016) andGitHub repositoryRatiometric control fortarget population (seeFigure S4D) α ( w opt /w T − α = 1, w opt = 1000 µ ; µ = 0 .
023 Control C in (Johnstonand Jones, 2016) andGitHub repositoryProduction independentof wild-type (see FigureS4E) α/w T ; α = 5 µ ; µ = 0 .
023 Control F in (Johnstonand Jones, 2016) andGitHub repositoryGeneral linear feedback(see Figure S4F) µ + b ( κ − δ w s − δ w f − δ m s − δ m f ); see TableS2, δ = 0 . δ = 1 . δ = 0 . δ = 0 . µ ; see Table S2 Control X in GitHubrepositoryRatiometric controlthrough degradation(see Figure S4G) λ ; λ = 0 . µw T /w opt ; µ = 0 . w opt = 200 Control G in (Johnstonand Jones, 2016) andGitHub repositoryLinear feedback in degra-dation (see Figure S4H) λ ; λ = 0 . µ + b ( w T + δm T − κ ); µ = 0 . b = 10 − , κ = 1000, δ = 1 Control Y in GitHubrepositoryDifferential control fortarget population indegradation (see FigureS4I) λ ; λ = 0 . α ( w T − w opt ); α = 1, w opt = 1000 Control Z in GitHubrepository38 igure S5. Ansatz predicts heteroplasmy variance for linear feedback control in a fast mtDNAturnover regime when fixation probability is low . Stochastic simulations of the linear feedback control networksystem with an mtDNA half-life of 2 days (Poovathingal et al. , 2012), corresponding to µ = ln(2) /
2. ( A ) The meancopy number, and ( B ) the mean heteroplasmy show qualitatively similar dynamics to the nominal parametrisationpresented (see Fig. 2C&D ). ( C ) Eq.(13) predicts V ( h ) accurately up to approximately 250 days, where the Eq.(13)begins to overestimate the variance. ( D ) The over-estimation of heteroplasmy variance coincides with an increase inthe probability of fixation at h = 0. Parameters apart from µ were chosen according to the protocol outlined in Choiceof nominal parametrization, with 10 repeats. κ = 101 . b = 2 . × − , for all other parameters see Table S2. x It e r a t e = wild-type= mutant= protected from death= susceptible to death Wait for time Death from suscetible population, birth from any individual
A B
Choose mutantsand wild-types to be susceptibleto death
Figure S6. Exploration of analogous Moran processes . ( A ) The original biallelic Moran process satisfiesEq. (S72), where h is the initial heteroplasmy, which is equivalent to E ( h ). ( B ) The “protected” Moran process. Thepopulation size is constrained to be fixed to some large constant, n . There exist two alleles, mutants (black circle) andwild-types (white circles). hnf s mutant and (1 − h ) nf s wild-type molecules are susceptible to death, the rest areprotected from death (denoted by a bar). An exponential random variable is drawn as the waiting time to the nextevent (see A modified Moran process may account for the alternative forms of heteroplasmy variance dynamics underdifferent models of genetic mtDNA control for discussion of the form of the rate Γ). Time is incremented by the waitingtime, then a death event occurs from the susceptible population and a birth event from any individual simultaneously.The same individual is allowed to be chosen for both birth and death. The process is then repeated iteratively. igure S7. A deterministic parameter sweep of fusion selectivity and the relative fusion rate formitochondrial quality control . An ODE treatment allows smaller heteroplasmy changes to be probed without theneed to resort to an infeasible number of stochastic simulations. Displaying the relative change in heteroplasmy (∆ h )after t = 1000 days. We observe that a reduction in heteroplasmy is achieved at intermediate fusion rates at allnon-zero fusion selectivities investigated. Grey denotes a change which is smaller than floating point precision. able S4. Comparison of previous models with our model of mitochondrial genetic/network dynamics .Each of the key differences with our model is enumerated, and has a corresponding comment; see Discussion. Model Key assumptions Key differences Comments (Tam et al. , 2013, 2015) • Cell consists of 16sub-compartments • Fission/fusioninduces migra-tion betweensubcompartments 1. Slower fission-fusion dynamicsresult in largerrate of increase in V ( h )2. Fast fission-fusionrates cause E ( h ) toincrease 1. Limiting toregimes wherefission-fusion isfast, likely re-sults in loss ofrate magnitudesensitivity2. Could be explainedby spatial effectswhich we neglecthere(Mouli et al. , 2009) • Fission follows fu-sion • Mitochondria con-sist of multipleunits ••