Energetic costs of cellular and therapeutic control of stochastic mitochondrial DNA populations

The dynamics of the cellular proportion of mutant mtDNA molecules is crucial for mitochondrial diseases. Cellular populations of mitochondria are under homeostatic control, but the details of the control mechanisms involved remain elusive. Here, we use stochastic modelling to derive general results for the impact of cellular control on mtDNA populations, the cost to the cell of different mtDNA states, and the optimisation of therapeutic control of mtDNA populations. This formalism yields a wealth of biological results, including that an increasing mtDNA variance can increase the energetic cost of maintaining a tissue, that intermediate levels of heteroplasmy can be more detrimental than homoplasmy even for a dysfunctional mutant, that heteroplasmy distribution (not mean alone) is crucial for the success of gene therapies, and that long-term rather than short intense gene therapies are more likely to beneficially impact mtDNA populations.

In this section, we use van Kampen's system size expansion (SSE) [1] to approximate the above equation. First, we derive the leading order solution of the expansion, known as the linear noise approximation (LNA). As we will see, the LNA provides useful insights and can accurately describe the dynamics of our system for short time-scales. However, for longer time-scales higher order terms are required to maintain accuracy.
Generally, a master equation can be written in the forṁ where Ω is the system volume, R is the number of reactions involved, N is the number of species, (S) ij is the stoichiometry matrix, n = (n 1 , n 2 , · · · , n N ) gives the number of particles of each species, E is a raising and lowering operator (its effect on an arbitrary function f ( n) is given by e.g. E i f ( n) = f ( n + e i ) and E −1 i f ( n) = f ( n − e i ) where e i is a column vector with all zeros and a 1 at entry i), and p n (t) is the probability distribution of states n at time t. The equations f j ( n, Ω) specify the reaction rates. The SSE provides an expansion of this master equation in powers of the inverse square root of the system volume Ω.
The expansion starts by expressing the state of the system in terms of a deterministic and stochastic component. Consider a single species whose concentration and copy number are denoted by φ(t) and n, respectively. One expects the distribution of n to be centered around Ωφ(t) and to have a width of order √ n ∝ √ Ω. This motivates the following equation where ξ i describe the fluctuations around the deterministic solution φ i (t). All terms in the master equation can now be expressed in the new fluctuation variables ξ i according to the following transformations: p n (t) = Π( ξ, t) Substituting Eqs. (4) into Eq. (2) gives which can be solved by collecting powers of Ω −1/2 . Terms proportional to Ω 1/2 form the macroscopic rate equations. The terms of order Ω 0 form a Fokker-Planck equation, the solution of which is the LNA: The coefficients A ij and B ij can be found by expanding Eq. (5) and are defined as The system we consider has two species (N = 2), four rate equations (R = 4) given by and a stoichiometry matrix by where the two rows refer to the wildtype and mutant species, respectively, and the columns represent the reactions f i . In order to continue, it is useful to rewrite Eq. (5) as an expansion in derivatives with respect to . . are functions of ξ, Ω and f .
2 Feedback control of a heteroplasmic mtDNA population

Steady states of a linear feedback control
When mutant and wildtype mtDNA molecules have identical replication (λ m = λ w ) and degradation (µ m = µ w ) rates, infinitely many steady states exist (i.e. states in which λ = µ). These steady states form lines in (w, m) space which can be straight (linear feedback control, Fig 1A), form segments of ellipses (quadratic feedback control) or take more complicated forms. Deterministic trajectories will asymptotically approach the steady state line while stochastic trajectories can fluctuate along the line thereby changing heteroplasmy. Similar lines of steady states are present in the relaxed replication model [3,4]. A linear feedback control λ(w, m) = c 0 − c 1 (w + δm) (with c 0 , c 1 > 0 constants) gives rise to a straight line of steady states, the slope of which is determined by δ. A small δ means that the steady state line intersects the mutant axis at higher copy number than the wildtype axis, meaning that total copy numbers are higher at h = 1 than at h = 0. In the extreme case of δ = 0, mutant copy numbers can fluctuate off to infinity, though in practice their numbers will be bounded by space restrictions and resource competition. When δ > 1 copy numbers will decrease as h increases; mutants are now sensed more than wildtypes, which could be caused by e.g. excessive production of reactive oxygen species by mutants (which is then sensed by the cell and incorporated in its feedback function). Fig. 1 in the main text shows the near-equivalence of mean wildtype, mutant, and heteroplasmy dynamics between three different forms of control, all of the form λ(w, m) = λ(w+δm): Fig. 1B, C, D shows that variance dynamics are also nearly identical.
We merely show that these different controls can give rise to similar dynamics by parameterising them such that they do. First, we set the variances of each control to be equal in the absence of mutants. From Table 1 it can be seen that the wildtype variance is now completely specified by the mtDNA degradation rate µ, the steady state copy number w ss and the control derivative ∂ w λ(w) evaluated at steady state. The latter quantity is given by i) −c 1 , ii) −d 1 /w 2 ss and iii) −2s 1 w ss for the three controls defined above. These expressions were all set to be −7.21 × 10 −6 , a value that was chosen such that the wildtype steady state distribution has a coefficient of variation of CV = 0.1. This value is arbitrary, other choices of give similar results. We use a deterministic wildtype steady state of 1000 (in the absence of mutants), which fixes the other control parameters. Note that it is not surprising that these different controls yield similar dynamics after we have parameterised them to do so; we want to illustrate these similarities to stress that which quantity is being controlled can be more important than how it is controlled. Other parameters used are Nss = 1000 (referring to the steady state copy number present in the absence of mutants), µ = 0.07 (corresponding to a half-life of 10 days), and initial copy numbers (w 0 , m 0 ) = (920, 160) (corresponding to an initial heteroplasmy of ∼ 0.15).

Homogenate heteroplasmy values
As mentioned in the main text, without explicit selection for either mtDNA species, mean cellular heteroplasmy remains constant at its original value, i.e. m0 m0+w0 where w 0 and m 0 denote the initial wildtype and mutant copy numbers, respectively. For long times, the fraction of cells that fixate on mutant species (i.e. h = 1) is given by m 0 /(w 0 + m 0 ) and these cells have mean copy number w opt /δ, whereas wildtype cells have mean copy number w opt . We can now calculate the homogenate heteroplasmy at long times (i.e. times at which all cells have fixated) as We note that Eq. 14 only holds when δ > 0.
3 The relationship between resource consumption and energy production: s(r i )

A linear output model
In our cost function we need to specify how s, the power supply measured in ATP/s (including leak) of a mitochondrion, depends on r i , a quantity resembling the resource consumption rate of a mitochondrion or type i (w or m). For robustness, we use two different equations for s(r i ) .
The first model is based on measurements in isolated mitochondria which found a linear relationship between r i and s (e.g. [5,6,7]), with r i referring to oxygen consumption rate. Specifically, we use the following equation: where φ can be mapped to the effective P/O ratio (explained in section 4.1) and β indicates the respiration rate at zero energy production and therefore specifies the amount of leak. We used the data from Ref. [5] to fit the parameters of this linear model ( Table 2,  It is not clear, however, whether the observations made in isolated mitochondria, without any interactions between the mitochondria and the nucleus or the endoplasmic reticulum, still hold in vivo. This is one of the reasons to also consider another type of model, as described below.

A saturating output model
Our cost function contains a term which penalises the consumption of resources, meaning that the cheapest state is the one that satisfies demand with the smallest net resource consumption rate. A linear function s(r i ) then implies that, for a given demand, the optimal number of mitochondria is the minimum number required to satisfy the demand with these mitochondria respiring as fast as they can (as can be seen in Fig 2E). However, containing a minimally required number of mitochondria will make a cell less robust to stochastic fluctuations in both mitochondrial copy numbers and demand. Moreover, mitochondria are known to have large spare capacities [8,9] indicating that in resting state they do not operate near their limits. We therefore expect that there is some extra cellular cost associated with this 'maximally respiring' state, causing it to be non-optimal in resting conditions. This is why we have chosen to use a second model which describes a saturating relationship between r i and s (visualised in Fig 2A). Note that we do not claim that mitochondria become less efficient as they respire faster, we impose the saturating shape merely to effectively assign a higher cost to high respiring states. By imposing this saturation, the variable on the y-axis of Fig 2A can be interpreted an 'effective energy production'. We will refer to the two models as 'the linear output model' and 'the saturating output model'.
A changing energy production efficiency is not entirely unreasonable, though, because in experiments with isolated mitochondria one usually uses a particular substrate (or a particular combination of substrates) whereas a larger mixture of substrates will be available in the cell, and the relative presence of each substrate may fluctuate over time. The efficiency of respiration depends on the kind of substrate that is used, so it may be possible that at high demand (and high respiration) the substrate of first choice has become limited and another less efficient substrate is used instead. Also, it was suggested that spare capacity can be caused by an increase in substrate entrance in the tricarboxylic acid (TCA) cycle [10]. This would mean that a state of high respiration is caused by an increase in electron transport chain substrates (e.g. NADH and FADH 2 ). A 'push' to the proton pumping complexes instead of a 'pull' at the ATP synthase would lead to an increase in the electrochemical gradient across the membrane and therefore an increase in leak. These arguments are speculative but show that the saturation model may not be unreasonable.
The saturating model is defined by the equation: where k is the rate at which the model saturates, s max is an indication of the maximum energy production rate, and ∆ > 0 is introduced to ensure that, like in the linear model, a nonzero amount of resource is required to maintain a zero energy output (due to proton leak). As proton leak accounts for ∼ 10% of the total energy production in the linear model, we chose ∆ = 0.1. We will set the parameter values k and s max such that the saturating and linear model show similar behaviour for low values of r w (section 4.5), as shown in Fig

Parameter values for the cost function
Some of our results will be a consequence of the exact structure of our cost function, and might have been different if another type of cost function was used. We would argue, however, that the main elements in our cost function are quite general: terms involving supply, demand, and resource. We aimed at making our cost function simple, and using biologically interpretable parameters. We do not aim to give a detailed kinetic description of the energetic costs involved, but present a simpler description that allows us to compare distinct strategies relative to each other rather than providing absolute costs. We use our cost function as a tool to characterise cost landscapes and begin to explore optimal control strategies.
In the spirit of 'back-of-the-envelope' reasoning in biology [11] we seek plausible and interpretable parameter estimates, using both order-of-magnitude estimations and values found in the literature. Default parameter values are summarised in Table 2.
The final goal is to obtain a cost for a state with a certain number of mtDNA molecules; we therefore need to express our cost as a cost 'per mtDNA molecule'. Because the density of mtDNA molecules within the mitochondrial network seems to be roughly constant [12], we assume every mtDNA molecule is associated to a particular amount of mitochondrial volume which we here define as a 'mitochondrial unit'. All our parameters refer to these mitochondrial units.

Mitochondrial energy production and leak: φ and β
Mitochondrial respiration is a process by which energy from nutrients is converted into (amongst others) ATP. Part of this process involves the pumping of protons across the inner mitochondrial membrane to create an electrochemical potential across the membrane. Energy is released when protons flow back into the matrix and this energy can be used to create ATP. However, the coupling between proton pumping and ATP synthesis is not perfect and protons can leak through the membrane, reducing the efficiency of respiration. An often measured quantity in experimental studies is the mechanistic P/O ratio [13,14] which refers to the theoretical maximum amount of ATP (P) produced per oxygen (O) reduced by the respiratory chain. The effective P/O ratio is more physiologically relevant and takes into account leak [7].
In our linear model we assume a linear relationship between mitochondrial ATP production rate and mitochondrial oxygen consumption rate, based on experimental data [5,7,6]. In Refs. [5,6] measurements show an almost perfect linear relationship between these two quantities, which is consistent with data from Ref. [7]. The parameters φ and β correspond to the slope and intercept of the linear function, respectively. In Ref. [5] this slope was measured to be 2.03 ± 0.13 for isolated pectoralis muscle cell mitochondria in the presence of pyruvate and malate; we have decided to use φ = 2, mainly based on these experiments.
The value of β is an indication of the 'leakiness' of the mitochondrion: it represents the rate of oxygen consumption that is required to balance the leakage of protons across the membrane in order to maintain the mitochondrial membrane potential. To obtain a consistent value for β we also use the data presented in Ref. [5]. Their measurements find that β is about a tenth of the maximum respiration rate. This maximum respiration rate is obtained by adding high (unlimited) concentrations of ADP; the state of the cell in these conditions is known as state 3 ADP . Because state 3 ADP does not necessarily correspond to in vivo conditions, we define the respiration rate in this state as r max,t : the maximum 'theoretical' respiration rate. We will fix r max,t = 1 and use this to scale our other parameters. This means that our parameter value for β is β = 0.1.
We stress that though we have based our parameter values here on a specific study, changes in their values will not affect the qualitative structure of our cost function but merely changes the slope and intercept of the linear output function defined in Eq. (15).

4.2
Resource and supply and maintenance cost: r max , r n and s n A cell in vivo is unlikely to experience the high concentrations of ADP present in state 3 ADP (defined above). We therefore set our parameter r max , the physiological maximum respiration rate, to be slightly below r max,t , i.e. r max = 0.95 · r max,t = 0.95.
We introduce r n and s n as the 'normal' respiration rate and ATP production rate which are present in resting conditions, respectively. Their values are not used directly in our model equations, but only to connect our parameter estimates to actual physiological values (section 4.8). Mitochondria have spare capacity, i.e. in normal unstressed conditions they use only part of their maximal oxygen consumption rate (OCR). The amount of spare capacity is usually measured as the fold-change in OCR that occurs after adding FCCP to cells, a mitochondrial uncoupler. Several measurements of the fold-change in OCR are: in the range (2-4)-fold [15], 1.4-to 2.5-fold [10], about 2-fold [8] and about 2.5-fold [9]. The maximum respiration rate when adding FCCP is known as State 3 FCCP , and is higher than State 3 ADP (see e.g. [16]). We interpret the spare capacity as the ratio State 3 FCCP /r n , meaning that the ratio r max /r n is lower than this. We have decided to take r max /r n = 1.5, meaning that r n = r max /1.5 ≈ 0.63. For the linear output model, the value of s n is now simply s n = s(r n ) ≈ 1.1 (this abstract value is related to actual values in ATP/s in section 4.8).

4.3
Maintaining, building and degrading: ρ 1 , ρ 2 and ρ 3 The model organism with the most quantitative data on mitochondrial energy budgets is budding yeast, Saccharomyces cerevisiae. Reasoning that the scales of biophysical costs of mitochondria are likely comparable across eukaryotes, we first draw from this literature to motivate order-of-magnitude estimates for a range of essential mitochondrial processes. We will later construct parallel estimates using other organisms.
We first focus on estimating the mitochondrial building cost ρ 2 (in ATP). We provide three distinct estimations and combine them to obtain our final estimate.
For the first estimation we use a list of mitochondrial proteins in yeast [17], and obtain information on turnover rates, abundance (per yeast cell), and lengths (amino acid length) of these proteins by using the Saccharomyces Genome Database [18]. We end up with a list of ∼ 200 mitochondrial proteins in yeast S. cerevisiae. We incorporate the observation that it takes about 5.2 ATP molecules to elongate a growing peptide chain by adding an amino acid [19], which means that the total synthesis cost of the mitochondrial proteins included in our list is given by 5.2 i length i abundance i ≈ 2 × 10 10 ATP per yeast cell. The known number of mitochondrial proteins in S. cerevisiae is on the order of 1000 [20]; we will therefore assume that the protein synthesis cost obtained from our protein list corresponds to roughly a fifth of the total mitochondrial protein synthesis cost. We might expect that the proteins best known [17] are the more abundant ones, meaning that our final cost is likely to be an overestimate. We also assume that all of the mitochondrially associated proteins are used exclusively for mitochondrial function. In E. coli, the mitochondrial protein synthesis cost represents ∼50% of the total mitochondrial synthesis cost (two other major contributors are phospholipid synthesis and RNA synthesis) [19]; we will make the assumption that this observation in E. coli holds in mitochondria as well. This brings the total mitochondrial building cost in a single yeast cell to be ∼ 1.9 × 10 11 ATP. Assuming 50-100 mtDNA molecules per yeast S. cerevisiae cell [21], the building cost associated with a single mtDNA molecule (and therefore the building cost of a mitochondrial unit) is given by (2 − 4) × 10 9 ATP.
For our second estimation we use the total protein weight of a single mitochondrion which was measured to be about 3 × 10 −10 mg in rat liver [22], as well as the typical weight of a single protein which is about 5×10 −17 mg [19]. This means that a mitochondrion contains about 6×10 6 proteins. Using that the typical length of a protein is 300 amino acids [23,24,19] together with the 5.2 ATP cost of adding amino acids and the estimation that protein costs represent 50% of the entire building cost, the mitochondrial building cost is estimated to be 2 × 10 10 ATP. Note that this cost does not necessarily represent our mitochondrial unit because it is unknown how many mtDNAs a 'mitochondrion' corresponded to when measuring its weight in Ref. [22].
The third estimation is based on the building cost of an E. coli, which is about 10 10 ATP [19]. Keeping in mind that we want the building cost of a fraction of mitochondrial volume corresponding to a single mtDNA molecule, we need to convert the building cost of an E. coli (which has a volume of about 1 µm 3 ) to represent our mitochondrial unit. It was estimated that the total mitochondrial network length in yeast S. cerevisiae is about 25 µm with a total mtDNA copy number of 50-100 [21]. Assuming that the mitochondria form tubules with a constant diameter of 300 nm [25] gives a total mitochondrial volume of about 1.8 µm 3 ; another total mitochondrial volume estimate in yeast S. cerevisiae is 1.5 µm 3 [25]. Uniform distributions for the mitochondrial volume (1.5-1.8 µm 3 ) and mtDNA copy numbers (50-100) leads to a volume of (2.3 ± 0.5) × 10 −2 µm 3 per mitochondrial unit. This means that rescaling the E. coli building cost gives us an estimate of (2.3 ± 0.5) × 10 −2 · 10 10 = (2.3 ± 0.5) × 10 8 ATP to build a single mitochondrial unit. Note that this may represent an overestimation because E. coli is a unicellular organism, whereas the mitochondrion is an organelle which cannot survive in isolation [26].
While there are differences in these estimates, arising both from uncertainty and different quantitative lines of reasoning, they together give an overall scale for mitochondrial building cost of around 10 9 ATP. Because in our model we only need a rough estimate of the mitochondrial building cost, we use ρ 2 = 10 9 ATP.
We interpret the maintenance cost, denoted by ρ 1 , as the cost in molecules ATP/s corresponding to, for example, maintaining the mitochondrial lipid membranes, importing/exporting proteins, and synthesizing new proteins. To obtain an estimation of ρ 1 , we again use the Saccharomyces Genome Database [18]. We can calculate the cost of continuously turning over the ∼ 200 proteins in our list (obtained from [17]), leading to 5.2 i length i abundance i (degradation rate) i ≈ 6 × 10 5 ATP/s per yeast cell. In other words, the maintenance cost per second of our set of mitochondrial proteins is about five orders of magnitude less than their synthesis cost. We therefore assume ρ 1 = 10 −5 ρ 2 .
The mitochondrial degradation cost ρ 3 is the most challenging parameter to estimate, as the process of mitochondrial degradation remains poorly characterised. Protein production and biosynthesis costs form the bulk of mitochondrial production requirements, and from cell-wide studies on energy budgets are among the most considerable demands in cell biology. We therefore assume that degradation has lower energy requirements than production, and set its upper limit at ρ 3 = 0.1 · ρ 2 . Lowering ρ 3 further has little impact on the outcomes of our model.

Demand and resource availability: D and R
We want to model two kinds of cells: low copy number cells (1000 wildtype mtDNAs in resting state) and high copy number cells (5000 wildtype mtDNAs in resting state). Denoting this desired number of 'normal' mitochondria by w n (which are assumed to operate at the 'normal' production rate s n (section 4.2)), we obtain where µ and λ are the degradation and replication rates per second. This equation states that the overall net output of w n mitochondria exactly satisfies demand. Using w n = 1000 and w n = 5000 leads to D ≈ 1055 and D ≈ 5275 ATP/s. These abstract values for demand are mapped to actual cellular ATP demands in section 4.8.
The parameter R denotes the maximum rate at which resource can be consumed by all of the mitochondria together and represents a cellular resource availability. In normal resting state the total resource that is consumed is w n r n , and the resource consumed when these mitochondria respire as fast as they can is w n r max . We then assume this maximal respiration rate is achieved by using all of the available resources, i.e. R ≈ w n r max . It may be, however, that a state of maximum respiration can only be maintained for a short time, and in our cost function we want to describe the 'steady state cost' for different states. We therefore use R < w n r max (R now denotes the maximal respiration rate that can be maintained for longer periods of time). We chose to use R = 0.8w n r max leading to R = 760 and R = 3800 for w n = 1000 and w n = 5000, respectively.

Saturating model parameters: s max and k
We set the parameters of the saturating model described in Eq. (16) such that it matches the linear model for low respiration. This lead to k = 3.0 and s max = 1.54.

The cost of resource consumption α
The value of α, i.e. the scaling parameter that appears in the cost function given by is hard to determine. Its value describes the cost of a unit of resource consumption relative to the cost of a unit of 'energy deficiency' (the cost of S(w, m) being one energy-unit below D). We estimated that a penalty for resource consumption usage should be about an order-of-magnitude less than the penalty for not-satisfying demand, and have therefore decided to assume α = 0.1. We note that the value of α has no influence on the shape of the demand-satisfying region, it only changes the relative costs within (and outside of) the region.

Mutant parameters: 1 , 2
In the main text we vary the parameter 1 , describing the resource uptake rate of mutants relative to wildtypes. Additionaly, mutants can be less efficient than wildtypes, producing less energy per resource consumed; we denote this lower mutant efficiency by 2 ∈ [0, 1]. Because the number of protons that are pumped across the mitochondrial inner membrane by the electron transport chain complexes for every unit of resource (NADH) that is consumed is fixed, a lower 2 would have to mean that either i) the mutation has increased proton leak (or other ways of depolarising the membrane), or ii) the mutation has made the ATP synthase dysfunctional. The value of 2 can be related to the P/O ratio of the mutants relative to that of the wildtypes. Most mtDNA mutations, however, affect the electron transport chain complexes themselves and are therefore likely to reduce the flow of resources through the chain (which would mean a low value for 1 ). This is why we assume 2 = 1 in our main model and only vary the parameter 0 < 1 < 1. For completion, here in the SI we provide a heatmap showing the cost in (w, m) space for various values of 2 (Fig 3).

Parameter units
We can relate our parameter values to actual values, e.g. expressing our demand D in ATP/s. As an example of an ATP demand we use the ATP production rate in unstressed human skin fibroblasts. Assuming that these healthy cells satisfy their demand, their net ATP production rate should equal their ATP demand. The rate of ATP production in skin fibroblasts was estimated to be about 10 9 ATP/s, the large majority of which is supplied by mitochondria [27]. The number of mtDNAs in healthy human skin fibroblasts was measured to be roughly in the range 2400-5200 [28] (the variation in copy number was partly due to variation in ages of the individuals), and we will use the value 4000 as an estimation.

2.0
β β denotes the resource consumption rate (in 'resource' per second) at zero energy production, and is therefore a measure of proton leak. Because we fix the maximum resource consumption rate to be 1 (section 4.1), β provides the fraction of resource consumption that is destined for 'futile' leakage. The mitochondrial building cost in units of 2.6 × 10 5 ATP. 3828 ρ 3 The mitochondrial degradation cost in units of 2.6 × 10 5 ATP. 383 1 The mutant energy production efficiency (the energy produced per unit of resource consumed) relative to that of wildtype mtDNA molecules. A low value of 1 can be caused by a high proton leak or a deficient ATP synthase.
The mutant resource uptake rate relative to that of wildtype mtDNA molecules. We will mainly use 1 = 1 and 2 < 1 because mtDNA mutations usually affect the proton pumping electron transport chain complexes. A defect in e.g. complex I will reduce its activity and therefore also the rate of consumption of NADH.   (16) for various parameter values. Note that none of the lines in the figure crosses the origin, because even when no ATP is created, respiration is required to maintain the gradient which would otherwise be lost due to proton leak. Fig 2B-E shows how resource consumption rate and cost change as the number of mtDNAs changes. Three regimes can be distinguished: 1) there are too few mitochondria present to satisfy demand and they use their maximum possible resource uptake to get their energy production rate as close to D as possible; 2) demand can be satisfied and the cost now only depends on the amount of resource that is used; 3) resource has become limiting and demand cannot be satisfied any more. The main difference between the two models is that in the linear model, when demand is satisfied, the cheapest state is the one with the smallest mtDNA copy number possible whereas the saturating model is cheapest at a higher copy number. This is because mitochondria in the saturating model becomes more efficient as less resources are being consumed per mitochondrion. Figs 3A,B show the cost function as a heatmap in (w, m) space. This figure is similar to Fig 2 in the main text except that here we have fixed 1 = 1 and varied 2 . Mutants are now less tolerated because they consume just as many resources as wildtypes, but still produce less output. It is now not the case that intermediate heteroplasmies are less efficient; intermediate heteroplasmies are only less efficient when 1 < 1 (mutants consume less resource) and when a saturating output model (Fig 2) is used. When 1 , 2 < 1, it is possible for intermediate heteroplasmies to be less efficient but the smaller the value of 2 , the smaller the range of values of 1 for which this is true; therefore, the effect of intermediate heteroplasmies being less efficient will be most easily observed when the mutant efficiency is close to that of the wildtypes. Figure 2: Relationship between resource consumption and energy output. A) The energy production rate of a single wildtype mitochondrion as a function of its resource consumption rate is shown, as given by Eqs. (15) and (16). For the linear model (corresponding to the straight lines) the parameters φ and β are changed by 10%, for the saturating model we vary smax and k. The magenta line indicates rmax. B) As w increases, demand is shared between more mitochondria and each individual one can afford to consume resources at a lower rate (the same figure legend applies for figures C, D and E). C) The total resource consumption increases with w because the mitochondria need to consume a non-zero amount of resources to produce a net energy output and each mitochondrion comes with a maintenance cost. D) The total energy produced by wildtypes increases when mutants are present. E) When demand is satisfied, the cost increases with w in the linear model, resulting in minimal costs when copy numbers attain the minimum number required to satisfy demand (1). In contrast, for the saturating model the cost decreases at first because as individual resource consumption drops, the energy production efficiency increases. Minimum cost now occurs when mitochondria are working most efficiently (2). Parameters 1 = 0.1 and 2 = 1.0 were used.

Expensive intermediate heteroplasmies -explanation and robustness
Here we attempt to provide an explanation for why, in certain cases, intermediate heteroplasmy values are more expensive than high or low values. In these cases, in high heteroplasmy regions it is more efficient to increase heteroplasmy even more, whereas in low h conditions it is more efficient to decrease h; this automatically implies there exists some intermediate heteroplasmy value that is least efficient. Figs 4A, B show the amount of resource consumed by individual wildtype and mutant mitochondria in four different example states. All states have identical total copy numbers (w + m = 10 4 ) and total outputs (equal to demand), but different heteroplasmy values (h = 0.1, 0.3, 0.7 and 0.9). When heteroplasmy increases, the individual resource consumption rates r w and r m both increase to compensate for the higher mutant copy number; this is true both in a low-h region (h increases from 0.1 to 0.3, Fig 4A) and a high-h region (h increases from 0.7 to 0.9, Fig 4B). However, the total resource consumption rate does not necessarily increase because the increase in h has caused a number of wildtype mitochondria to become mutants, thereby decreasing the combined wildtype resource usage. Computing the values of the resource consumption rates for the states with h = 0.1 and 0.3, while referring to Fig 4A, gives: w 1 r w1 = 9000 · r w1 ≈ 3227 m 1 r m1 = 1000 · r m1 ≈ 125 w 2 r w2 = 7000 · r w2 ≈ 2981 while the states h = 0.7 and 0.9 (referring to Fig 4B) give Comparing the states with h = 0.1 and h = 0.3, the total rates of resource usage are 3352 and 3428 respectively; the lower heteroplasmy state is more efficient. However, when heteroplasmies are higher, the high heteroplasmy state (0.9 rather than 0.7, with total resource usage 3511 vs 3600) is most efficient. This effect is due to the nonlinearity of Eq. 16. Next, we investigate how robust the existence of an expensive intermediate heteroplasmy value is. Let h max denote the value of h with maximum cost at fixed total copy number (w + m). Fig 5A shows h max as a function of both total copy number and 1 , setting all other parameters to their default values. The range of total copy numbers was chosen to correspond roughly to the range for which demand can be satisfied in a fully wildtype cell (with the objective of considering 'plausible' copy number regions). We see that   B) The parameter smax is increased by 50% with minimal effect on the output. We kept the parameter k fixed as it defines the amount of proton leak (the resource consumption rate at zero energy output) which agrees with the amount of proton leak in the linear output model whose parameters are based on experimental data. C) The parameter ρ 1 is increased by an order of magnitude with minimal effect on the output. D) The parameter ρ 1 is increased by an order of magnitude and smax is decreased by 50%. Again, change in hmax are small.

Comparison of the cost of different control mechanisms
In the main text we introduced four different feedback controls and compared their mean costs. Here, we additionally show the means and variances of the wildtypes, mutants and cost up to ∼ 82 years resulting from stochastic simulations. The parameter values we used to compare the four controls are summarised in Table 3 for clarity. We observe that the increase in cost for the linear control with δ = 0 is mainly caused by an increase in mean mutant copy number.  Table 3. Again, we see that the effects of the control are more noticeable in low copy number cells. Parameter are set as given in Table 3. Values for wopt are those for the saturating output model at low and high copy number. The free parameters in control III and IV (δ and η) were optimised over initial conditions in the range h ∈ [0, 0.2]. For the optimization the default cost function parameters were used as well as 1 = 0.3. Fig 7 shows the optimal values for δ in the linear feedback control λ(w, m) = µ + c 1 (N ss − (w + δm)) as a function of 1 . Stochastic simulations starting in steady state at either h 0 = 0.1 or h 0 = 0.8 were performed for T = 10 4 days. The mean integrated cost over these 10 4 days was evaluated for different values of δ, and the optimal δ values are shown. This was done for both the linear and the saturating model. The general trend is, as was shown for T = 100 in the main text, the lower 1 the lower the optimal δ. Label Control optimal parameters satisfying our two constraints Table 3: Parameter values for the four different control mechanisms we employ. Two parameters of each control are set by the two constraints we impose. The parameter α R was proposed to lie in the range 5-17 [3] and here we used α R = 10. The values for δ and η are found by optimizing our cost function over the steady states corresponding to our initial conditions. We used 50 initial conditions equally spread over the range h 0 ∈ [0, 0.2]. The two values used for wopt are 1524 and 7616 (Table 2). We further use µ = 0.07 day −1 .

Heteroplasmy values can increase after nuclease treatments
As mentioned in the main text, there is a possibility for cellular heteroplasmy to increase after a treatment has been applied. This is true especially if the selectivity of the treatment is low (i.e. ξ is close to 1) and the initial heteroplasmy of a cell is high; in this case treating a cell may even eliminate all wildtype mitochondria, increasing heteroplasmy to 1. To model the extent of this effect we initialise a cell with a given heteroplasmy h 0 , and let it undergo one round of treatment and recovery after which the final heteroplasmy is recorded. This process is repeated to obtain the probability that, given an initial heteroplasmy h 0 , the final heteroplasmy after treatment exceeds h 0 (P (h final > h 0 |h 0 )). Figs 8B,C show these probabilities as a function of h 0 and selectivity parameter ξ, for initial mtDNA copy numbers 500 and 5000. We observe that the effect-size is larger for low copy number cells. Fig 8D shows an example of the distribution of post-treatment heteroplasmies. The recovery time used in the simulations is 30 days which is long enough for the cells to recover their initial copy numbers and short enough for the change in h to be almost completely due to treatment, rather than due to naturally occurring random drifts in heteroplasmy values. Chances of increasing heteroplasmy are highest when h 0 is very high or very low (if h 0 is low the low mutant copy numbers increase the effect of stochastic fluctuations).
In the examples shown, when ξ <= 0.6 (i.e. for every mutant that is cleaved, 0.6 wildtypes are cleaved) increases in heteroplasmy are very unlikely to occur.

Measurement of mtZFN expression profile
It was previously shown that the mtZFN pairing NARPd(+) and COMPa(-) specifically cleaves mtDNA with the m.8993T>G mutation [29,30,31]. We confirmed the transient expression profile of mtZFNs used in the modelling by transfecting 143B cells with a plasmid encoding NARPd(+) that co-expresses fluorescent marker protein mCherry. At 24 hours post-transfection, transfected cells were sorted using fluorescence activated cell sorting (FACS), to ensure a homogenously transfected sample (Fig 9A). Cells were harvested at 24 hours, and the remainder returned to culture dishes, to be harvested at later time points. Total protein was extracted from cells and analysed by western blotting, using antibodies to the HA epitope. Coomassie (CM) staining was used to verify equal loading of the gel. For full protocol see Ref. [32]. We observe that expression of mtZFNs is almost undetectable by 96 hours, and is totally undetectable by 120 hours (Fig 9B).

MtDNA copy number measurements in heteroplasmic cells
To determine the absolute mtDNA copy number in the pre-treatment 80% heteroplasmy cells, a plasmid construct containing the mitochondrial region of interest was created. The exact molecular weight of the dsDNA plasmid molecule was calculated based on nucleotide composition (3318.8 kDa). The plasmid and experimental sample concentrations were measured in triplicate with a Qubit fluorometer using the Qubit dsDNA BR Assay Kit (ThermoFisher Scientific). Three 10-fold plasmid serial dilutions were made (spanning 8 orders of magnitude) and analysed by quantitative PCR. Five similar dilutions of the experimental sample were measured in triplicate. Highest and lowest values were removed from further analysis to reduce the effect of pipetting errors. The calculated mtDNA copy number per µL sample (1.13 × 10 7 ) was determined as the average number for each unknown replicate using the average result based on each of the three standard curves. The number of cells per µL was calculated based on the mass of DNA per single diploid cell (6.57 pg) [33]. As 143B cells are mildly aneuploid, this is likely to be a small underestimate.
Primers and probe used for qPCR: mtDNA 3211 F: 5 CACCCAAGAACAGGGTTTGT 3, mtDNA 3298 R: 5 TGGCCATGGGTATGTTGTTAA 3 mtDNA 3242 probe: 6-Fam-TTACCGGGCTCTGCCATCT-Tamra Using this method, we calculate that each pre-treatment 143B cell (with heteroplasmy 0.8) in the initial cycle of iterative mtZFN treatments contained 889 ± 214 (S.E.) mtDNA copies per cell. Based on these measurements, we use a total initial copy number of 900 in our MCMC simulations used to fit the experimental data provided in Ref. [31].

Bayesian inference model outputs
Therefore, as long as the ratio I 0 /b is constant, a larger value for b does not influence the mtZFN dynamics and therefore the mtDNA dynamics. Because it is the ratio I 0 /b that determines the mtZFN dynamics at large b, we have performed the inference using this ratio rather than the parameter I 0 itself. Fig 11 shows that our model predictions of mtZFN expression profiles are in broad accord with the experimental data obtained in this study. Our prior distribution on the treatment duration parameter b allows for much more slowly decaying mtZFN concentrations, but our inference selected values for b that predict very low mtZFN concentrations 5 days post-transfection, agreeing with the experimental data.
Prior distributions used are: δ: N (1.0, 1.0) bounded between 0 and 10; log 10 I 0 /b: U (−1, 3); log 10 b: U (−3, 3); ξ: U (0, 1); log 10 c 1 : U (−7, −2); σ h : U (0, 0.5); σ T : U (0, 0.5). The priors of the standard deviations were chosen to be U (0, 0.5) because both h and T mostly lie in the range [0, 1] and experimental data indicate a standard deviation smaller than 0.5. The prior distribution of c 1 was motivated by investigating the general behaviour of mtDNA dynamics for various values of c 1 . We found that if c 1 = 10 −2 the control is unrealistically strong whereas if c 1 = 10 −7 it takes unrealistically long time to return to a steady state value when out of equilibrium. Fig 6 in the main text Figs 6A,B in the main text were obtained using stochastic Gillespie simulations in which mtDNA dynamics were simulated using Eq. (12) in the main text. 10 4 cells were simulated for 400 days with the following initial heteroplasmy distributions, all with mean h = 0.8: i) all h values fixed at 0.8 (left), ii) Beta(12, 3) (middle), and iii) Beta(0.4, 0.1) (right). The following parameters were used: w opt = 5000, I 0 = 39.6, b = 12.4, ξ = 0.76, c 1 = 5.1 × 10 −5 , δ = 1.0 and µ = 0.07 day −1 . These parameter values were chosen because they provide good fits to the experimental data in Ref. [31]. Fig 6C is obtained from stochastic simulations using our default cost function parameters for 'low' copy number settings (Table 2), and using the following mtZFN treatment parameters (which provided a good fit to experimental data): b = 11.9, ξ = 0.76, c 1 = 2.5 × 10 −4 , δ = 1.0 and µ = 0.07 day −1 . For each value of 1 , the value for I 0 corresponding to the minimum cost over a simulation time of 400 days was found.

Parameter values used in
The short and strong treatment in Fig 6D (blue line) uses b = 11.9, ξ = 0.76, c 1 = 2.5 × 10 −4 , δ = 1.0, and µ = 0.07 day −1 . We further used I 0 = 47.1 which was found to be optimal for 1 = 0.2 (the cost function corresponding to this value of 1 is shown in the figure). The long and weak treatment uses b = 0.1 and I 0 = 0.6 (the optimal treatment strength when using 1 = 0.2 and b = 0.1). For the more selective treatment (magenta line) we set ξ = 0.4, b = 11.9 and I 0 = 41.75 (the optimal treatment strength under these ξ and b). The heteroplasmy mappings corresponding to short and long treatments shown in Fig 6E were obtained from stochastic simulations using parameters identical to the ones for the short and long treatments in Fig 6D. We note that in finding I 0,opt , we initialised all cells at identical states. When a distributions of initial states is used, the variance that is now present is likely to affect the optimal treatment strength.   Fig 9B and were subsequently rescaled to investigate whether our model can broadly account for the experimentally observed dynamics (our predicted mtZFN concentrations are proportional to the measurement data points with an arbitrary proportionality constant).