## Figures

## Abstract

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.

## Author summary

Mitochondria, best known for their role in energy production, are crucial to the survival of most of our cells. To respond to energetic demands and mitigate against mutational damage, cells control the mitochondrial populations within them. However, the character of these control mechanisms remains open. As experimental elucidation of these mechanisms is challenging, theoretical approaches can help us understand the general principles of cellular control of mitochondria in physiology and disease. Here, we use stochastic modelling to compare control strategies by studying their impact on the dynamics of mitochondrial DNA (mtDNA) populations as well as their energetic burden to the cell. We identify optimal strategies for the cell to control against mtDNA damage and preserve energy production and use this theory to explore the action of recently developed mitochondrial gene therapies, which reduce the fraction of mutant mtDNA molecules inside cells. We show how treatment efficiency may depend on pre-treatment distributions of mutant and wildtype mtDNA molecules: treatments are less effective for tissues consisting of cells with highly varying mutant levels, and long-term, rather than short intense, gene therapies should be favoured.

**Citation: **Hoitzing H, Gammage PA, Haute LV, Minczuk M, Johnston IG, Jones NS (2019) Energetic costs of cellular and therapeutic control of stochastic mitochondrial DNA populations. PLoS Comput Biol 15(6):
e1007023.
https://doi.org/10.1371/journal.pcbi.1007023

**Editor: **Daniel A. Beard,
University of Michigan, UNITED STATES

**Received: **October 13, 2018; **Accepted: **April 11, 2019; **Published: ** June 26, 2019

**Copyright: ** © 2019 Hoitzing et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability: **All relevant data are within the manuscript and its Supporting Information files.

**Funding: **The author(s) received no specific funding for this work.

**Competing interests: ** The authors have declared that no competing interests exist.

## Introduction

Most human cells contain 100-10,000 copies of mitochondrial DNA (mtDNA) which are situated inside the mitochondria. The proteins encoded by mtDNA are crucial for mitochondrial functionality, and mutations in mtDNA can cause devastating diseases [1–6]. Heteroplasmy, the proportion of mutant mtDNA molecules in a cell, typically has to pass a certain threshold (∼ 60-95%) before any biochemical defects can be observed [7–14]. The existence of thresholds at which mutant loads begin to have an effect has profound implications for our understanding of disease onset, drawing attention to the variance dynamics of the mutant fraction in cellular populations. As this variance increases more cells can be above threshold, and thus show pathology, even if average mutant load is unchanged.

Mitochondrial biogenesis and maintenance require cellular resources, and mitochondria are key sources of ATP and play other important metabolic roles. The particular ‘effective cost’ that cellular control of mitochondria acts to minimise remains poorly understood: for example, both decreases [15] and increases [15, 16] in wildtype copy numbers have been observed for different mutations as the mutant load increases. Some studies suggest that mtDNA density is controlled [17–19], others that total mtDNA mass [20, 21], or mtDNA transcription rate [22] is controlled. Understanding mtDNA population dynamics inside cells, and how these populations react to clinical interventions, is crucial in understanding diseases [23, 24]. However, experimental tracking of mtDNA populations over time is challenging, necessitating predictive mathematical modelling to provide a quantitative understanding.

In parallel with efforts to elucidate cell physiological control, protein engineering methods to artificially control mtDNA heteroplasmy are making fast progress. Two recently developed methods for cleaving DNA at specific sites involve zinc finger nucleases (ZFNs) and transcription activator-like effector nucleases (TALENs) [25–31], which have been re-engineered to specifically cleave mutant mtDNA [32–36]. MitoTALENs have been successfully used to reduce mutant loads in cells containing disease-related mutations, but elimination of the target mutant mtDNA was not complete [32, 37]. Similarly, treating cells multiple times with mtZFNs led to near-complete elimination of mutant mtDNAs [35, 36]. Quantitative theory for these therapeutic technologies has not yet been developed, leaving open questions about how these tools can be optimally deployed.

In this paper, we develop theory from bottom-up bioenergetic principles which allows us to study the effects of distinct cellular mtDNA control strategies, to analyse the bioenergetic cost of different mtDNA states, and to combine mtDNA control and energy-based cost to identify optimal control strategies for the cell. Finally, we construct a model for therapeutic mtDNA control using recent experimental data [36] and highlight challenges linked to heteroplasmy variance.

## Results

### Control: General insights on the role of feedback control

We employ a linear form of mtDNA feedback control and assume each mtDNA molecule replicates and degrades according to Poisson processes with rates λ and *μ*, respectively. Because control of biogenesis or autophagy yield similar behaviours [38], we assume that the degradation rate *μ* is constant and that feedback control is manifest through the replication rate λ(*w*, *m*), where *w* and *m* denote the number of mutant and wildtype mtDNA molecules in the cell. To connect with experiments, we use *μ* ≈ 0.07 day^{−1} corresponding to a half-life of about 10 days [39]. We only model post-mitotic cells, though our analysis can be extended to include cell divisions.

Specifically, we use a birth rate of the form:
(1)
where *c*_{1} > 0, *w*_{opt} > 0 and *δ* are constants, with *w*_{opt} denoting the steady state value towards which the effective population, here defined as *w* + *δm*, is controlled. The magnitude of *c*_{1} determines how tightly the population is controlled. We use the term ‘mitochondrial sensing’ to describe how the cell might sense the mitochondrial population that is present. ‘Mutant sensing’ then refers to how strongly mutants are sensed relatively to wildtypes, which is encoded in the parameter *δ*. When steady state is reached (i.e. *w* + *δm* = *w*_{opt}), replication and degradation rates are equal. In the absence of mutants, the resulting wildtype steady state is assumed to be optimal. We note that assuming the existence of *w*_{opt} does not imply a control based on copy number. Other quantities related to mitochondria may be controlled instead, such as total mitochondrial mass or ATP production, their desired values being reached at an effective population size of *w*_{opt}. Thus, we define ‘mitochondrial sensing’ to refer to a wide range of mechanisms available to the cell to infer properties of its mitochondrial population, which can then be used to decide on a control action.

The deterministic dynamics resulting from this control are described in Eq (4). We do not include the possibility of *de novo* mutations but our approach can straightforwardly describe the subsequent behaviour if new mutations arise. Our linear model shares features with the ‘relaxed replication model’ [40, 41] (Eq (5)), though is written in a simpler form. The relaxed replication model has been used in a variety of other models [42, 43] and has obtained experimental support [15].

We will first investigate properties of more general control strategies, after which we return to our linear control and discuss parameterisations that optimise the energy status of the cell. Finally, we use the linear control to fit recent experimental data involving treatment of heteroplasmic cells with mtZFNs.

#### A wide range of control strategies induces similar mtDNA behaviour and admits quantitative analysis.

Many possible control strategies can be parameterised to give rise to nearly identical means and variances for wildtype, mutant, and heteroplasmy dynamics up to long times (e.g. a human life-time) (Fig 1A, 1B and 1C, S1B, S1C and S1D Fig). This is especially true when mtDNA copy numbers fluctuate around their steady state values, in which case a linear control forms a good first order approximation to the complex ‘true’ control function. In this case, it is not the manner in which the controlled quantity is being controlled, but **which** quantity is controlled that is the most important. For example, the extent to which mutants and wildtypes contribute to the replication feedback function (determined by *δ*) determines how their relative means and variances evolve (Fig 1D), and contains more information on the dynamics than the functional form of λ(*w*, *m*) for fixed *δ* (e.g. whether λ(*w*, *m*) is linear or quadratic).

**A-C: A wide range of cellular control strategies can yield similar dynamics**. Stochastic simulations were used to compare three structurally distinct cellular controls (see legend), each reflecting a different function of the underlying sensed quantity *w* + *δm* with *δ* = 0.5. All controls are set to have the same wildtype mean and variance in the absence of mutants (section 2 in S1 File). No explicit selection for either mtDNA species is used. Figure (C) illustrates the difference between cellular mean and tissue homogenate heteroplasmy. **D: Control tradeoffs are required when multiple species are present**. The more strongly one species is controlled, the more control is lost over the other. Changes in variances between cells as described by the linear noise approximation (section 1 in S1 File) are shown (intermediate times). For long times, fixation occurs and the variance of the surviving species saturates. In the most-right figure we depicted the case in which mutants contribute less to the control than wildtypes (*δ* < 1).

We stress the difference between two types of average heteroplasmy, as was also stressed in Ref. [41]: the individual cellular mean heteroplasmy (with *n*_{cells} the number of cells in the tissue) and the tissue homogenate heteroplasmy . This difference is clearly seen in Fig 1C. When no explicit selection is present for either mtDNA species, mean cellular heteroplasmy remains constant at its initial value *m*_{0}/(*m*_{0} + *w*_{0}), where *w*_{0} and *m*_{0} denote the initial wildtype and mutant copy numbers, respectively. The homogenate heteroplasmy at long times is given by *m*_{0}/(*m*_{0} + *δw*_{0}) (section 2.2 in S1 File). It is clear that when mutants contribute little to the feedback control (small *δ*), tissue homogenate heteroplasmy can reach high values, and even approach 〈*h*〉_{homog} = 1, without explicit selection. A tissue can thus appear, when studying the homogenate heteroplasmy, to show selection for one type of mtDNA over another, whereas in fact mean cellular heteroplasmy is unaltered and mutant and wildtypes have identical proliferation rates.

#### Nonlinear cost functions predict changes in tissue maintenance.

The birth-death model used to describe mtDNA dynamics can be written as a master equation (Methods, section 1 in S1 File). S1 Table shows the first order solution of the system size expansion, an approximation method to master equations, which is known as the linear noise approximation (LNA). Applying this approximation to a general form of mtDNA control (section 1 in S1 File), we find that i) if only one species is controlled, the variance of this species quickly reaches a constant value (see also [38]), ii) when both species are controlled with equal strength their variances increase at identical rates, iii) in general the more tightly controlled species has a more slowly increasing variance, and iv) the rate of increase of heteroplasmy variance depends, to first order, only on mtDNA copy number and turnover (as found in [38]) (Fig 1, Table 1()). Eventually, all variances reach a constant value due to fixation.

What are the biological implications of these findings? A given mtDNA state (*w*, *m*) will accrue a cost to the cell, denoted by *C*(*w*, *m*), which can e.g. be an energetic cost or some other metric of tissue burden. If this cost function is nonlinear, increasing variances in *w* and *m* can lead to changes in mean cost 〈*C*(*w*, *m*)〉 even when mean cellular copy numbers 〈*w*〉, 〈*m*〉 remain constant (Methods) because the mean of a nonlinear function of random variables is not generally equal to the function of the mean of those variables (as seen above with cellular vs homogenate heteroplasmy). Therefore, the mean cost of maintaining a tissue may increase over time, even if tissue demands and mean mtDNA levels stay constant (Table 1()). However, these increases may be small and their significance depends on the details of the cost function: hence the need to consider explicit forms, as we do in the next section.

### Cost: An effective mitochondrial energy-based cost function

Next, to find general quantitative principles underlying mitochondrial energy budgets, we build a cost function that assigns a cost to any given mtDNA state (*w*, *m*) and allows a general quantitative investigation of the tradeoffs in maintaining cellular mtDNA populations. The ‘true’ energy budget of a cell with a given mitochondrial population is highly complex, involving many different metabolic processes in which mitochondria are involved [44–46]. We provide a simpler description, focussing on ATP production as a central mitochondrial function, and removing kinetic details in favour of a coarse-grained representation, to provide qualitative rather than quantitative results.

#### General cost function structure.

Three important terms involved in the energy status of a cell are: i) the energy demand *D*, ii) the net energy supply *S*, and iii) the efficiency of energy supply. Here we define efficiency as the amount of energy produced per unit of resource consumed. We included intuitive and general terms in our energy-based cost function, such as replication, degradation and maintenance costs, supply and demand, and resource availability. We seek a cost function that captures the idea that there might be an optimal number of mitochondria: not so few that each mitochondrion is inefficiently overworked and not so many that the burden of the mitochondrial population is itself large.

We express our effective cost function as:
(2)
where *α* is a scaling constant, and *r*_{i} gives the rate of resource consumption of a mitochondrion of type *i* (*w* or *m*). The second term assigns a cost to the use of resource. The terms in this cost function are expressed as rates: *S* and *D* correspond to net energy production (supply) and demand per unit time. Supply and demand terms are left deliberately generalisable to encompass the differences in metabolic poise between cell types. The demand can be considered to represent energy requirements of all cellular processes besides mitochondria (whose maintenance costs are incorporated in their net supply *S*(*w*, *m*)), which we assume to be constant. We are therefore modelling post-mitotic cells in stable environments, as demands are expected to change throughout the cell cycle. This cost function can be evaluated for any state (*w*, *m*) and assigns the lowest cost to a state that satisfies demand in the most efficient way.

The net energy production of a state (*w*, *m*), *S*(*w*, *m*), is modelled as
(3)
where *ρ*_{1,2,3} are mitochondrial maintenance, building, and degradation costs, *s*(*r*_{w}) denotes the power production (in ATP/s) of a single wildtype mitochondrion given a resource consumption rate *r*_{w} (Methods, section 4 in S1 File), and λ and *μ* (which can be functions of *w* and *m*) denote the birth and death rates in units per second. Mutant mtDNA molecules are distinguished by the parameters *ϵ*_{1}, *ϵ*_{2} ∈ [0, 1] describing the mutant resource uptake rate (*ϵ*_{1}) and efficiency (*ϵ*_{2}) relative to that of the wildtypes (*r*_{m} = *ϵ*_{1}*r*_{w}). A low *ϵ*_{1} could represent reduced flow through the electron transport chain due to e.g. damaged respiratory complexes, whereas a low *ϵ*_{2} could denote increased proton leakage. Many mutants are known to have dysfunctional respiratory chain complexes [47], likely causing reduced electron flow through the respiratory chain and therefore reduced consumption rates of respiratory substrates such as NADH and oxygen. We therefore use *ϵ*_{1} < 1 and *ϵ*_{2} = 1 as our default choice (further described in section 4.7 in S1 File), though settings with *ϵ*_{2} < 1 are also investigated (section 5 in S1 File).

The more detailed structure of our cost function, which has been relatively general so far, comes from specifying the relation *s*(*r*_{w}). In other words, how does the energy output of a mitochondrion depend on its resource (e.g. oxygen) consumption rate? We consider two possible forms for this function (section 3 in S1 File): a linear output relationship, suggested by some literature [48–50], and a saturating relationship, which accounts for finite resource consumption and spare mitochondrial capacity [51, 52] (section 5 in S1 File). We refer to these alternatives as the ‘linear output model’ and the ‘saturating output model’. Both models are described in more detail in our Methods section (Eq 7).

A trade-off between yield (efficiency) and rate of ATP production is present in yeast [53–55] whose rate of ATP production due to respiration can become saturated at high resource levels or limited oxygen supply [53, 56]. Higher energy production rates can then still be obtained by using fermentation at the expense of a lower yield [53]. A similar trade-off may exist in mitochondria, whose power production efficiency is higher when oxidizing NADH compared to oxidizing succinate [50]. The former may be the preferable substrate due to its higher yield, but if its levels become limiting an increase in the relative use of succinate would lower overall efficiency. When oxygen is limiting, increased glycolysis in an attempt to increase ATP production also leads to lower overall efficiency. These mechanisms could be the cause of a reduced power production efficiency at high energy demand, as proposed in our saturating output model. We will contrast our findings of the saturating output model with those generated by the linear output model.

We fitted the parameters of the linear output model using data provided in Ref. [48] (section 3 in S1 File), and set the parameters of the saturating model such that the two models behave similarly at low *r*_{i}. Further details on the choice of parameter values, and their biochemical interpretations, are given in section 4 in S1 File; default values are provided in S2 Table. For a given demand *D* we find the *r*_{w} which gives a demand-matching supply (details of the scenarios when supply cannot meet demand are given in section Methods). Given *r*_{w} we calculate *r*_{m} (using *r*_{m} = *ϵ*_{1}*r*_{w}) and so calculate the cost *C*(*w*, *m*).

In taking both the linear and saturating output models into consideration we have endeavoured to build the most general picture of a mitochondrial cost function that retains bottom-up interpretability. Where possible, we estimate associated parameter values based on experimental data. However, other cost function choices are certainly possible and can be analysed using the platform we present below: our objective here is to complement the generic result regarding cost functions in paragraph “Nonlinear cost functions predict changes in tissue maintenance.” with a specific reasonable choice of cost.

#### Intermediate heteroplasmies may be inefficient and resource availability can dictate the cost of mtDNA states.

Fig 2 shows heatmaps of the cost function in (*w*, *m*)-space for different mutant pathologies (modelled as different values of *ϵ*_{1}). Our cost function generates a heteroplasmy threshold, its value depending on both *w*_{opt} and *δ*, above which demand cannot be satisfied using oxidative phosphorylation (Fig 2), though increased glycolysis may still maintain cell viability. The threshold effect is an established phenomenon in mitochondrial physiology [7–14].

A visualization of the cost function in (*w*, *m*) space is shown for both saturating and linear output models, for various mutant pathologies (described by *ϵ*_{1}). For visualization purposes, states in which cellular demand cannot be satisified are shown in white. Cells in these states may still survive by e.g. increasing glycolysis (effectively reducing mitochondrial demand). This figure assumes high copy numbers, results are qualitatively similar for low copy numbers. The actual cost values (given by the colour map) are of lesser importance for our findings, we rather focus on the qualitative shape of the cost function. **A**: The magenta (solid) and black (dashed) lines show the contour of the demand-satisfying region when demand is increased by 10%, or demand is increased by 50% and cellular resource availability is increased by 35%, respectively. **B**: The orange line corresponds to constant total copy number; moving up along this line increases heteroplasmy. Cells in region 1 or region 3 are more efficient, and show a lower cost, than cells in region 2. **C**: The linear mitochondrial output model does not show a decreased efficiency at intermediate heteroplasmy values.

The state with lowest cost according to the linear output model is one with a minimum number of mitochondria required to satisfy demand (Fig 2C), where these mitochondria respire as fast as possible. This would mean that this state of lowest cost has no spare capacity. Assuming a cell controls its mitochondrial population towards the state with lowest cost, the linear model predicts cells to lack spare capacity, contradicting experimental observations [51, 57, 58]. A saturating output model solves this problem (S2E Fig and Fig 2A and 2B) and generates a trade-off between using each mitochondrion efficiently (minimising its resource consumption by increasing population) and minimising the cost of maintaining the total number of mitochondria (achieved by reducing the population). At low resource consumption, representing the linear regime of the saturating output model, the two models are similar.

We observe other qualitative differences in cost function structure between the saturating and linear output models. In the former, it is possible for intermediate heteroplasmy states to be more expensive than states homoplasmic for either species (Fig 2A and 2B). Hence, in the saturating output model, it is possible for intermediate heteroplasmies to be the least efficient and the most expensive (Fig 2 and S4 Fig, Table 1()). This result arises from a tradeoff, when mutant load is increased, between a decrease in global efficiency and a reduction in resource consumption by the new mutants (section 5 in S1 File). The linear output model, on the other hand, always shows higher costs at higher heteroplasmy (for fixed total copy number) (Fig 2C).

For our cost function, the existence of a high-cost intermediate heteroplasmy value is a relatively general feature of the saturating output model. We calculated the value of heteroplasmy with maximum cost (at constant total copy number), denoted by *h*_{max}, as a function of several model parameters. *h*_{max} = 1 (the homoplasmic mutant state) at values *ϵ*_{1} ≲ 0.3 due to very low mutant functionality. However, at higher values of *ϵ*_{1} (0.5 ≲ *ϵ*_{1} < 1) we find *h*_{max} ∼ (0.5–0.8) over a large range of several of our cost function parameters (section 5.2 in S1 File). Though the size of the effect may be small, its existence alone is an interesting feature of our saturating output model.

It was previously found that it is possible for two mtDNA variants in mice to function normally at homoplasmy, but show deficiencies in heteroplasmic states [59]. While we do not claim that our model is the reason behind these observations it does suggest that differing resource consumption rates associated with distinct mtDNA species may play an important role.

### Combining cost and control: Comparison and optimisation of both cellular control and treatment strategies

#### Timescales and energy sensing in optimal control of mtDNA populations.

Here we compare the mean cost over time for four plausible cellular control strategies. The first two consist of the linear feedback model λ(*w*, *m*) = *μ* + *c*_{1}(*w*_{opt} − (*w* + *δm*)) with (I) *δ* = 0 (only wildtypes are sensed) and (II) *δ* = 1 (total mtDNA copy number is controlled). We further identify optimal parameterisations (i.e. ones that minimise steady-state cost) of two control strategies, namely (III) a linear feedback control and (IV) the ‘relaxed replication model’ (Eq (5)) [40, 41].

First, we fix the parameter values that are not being optimised. Our goal is to compare the costs of the dynamics resulting from each of the four controls, in the presence of mutants. In other words, we want to investigate whether some of these controls are better at protecting the cell (in the sense of maintaining a low energy cost) against mutant loads than others. To make this comparison fair, we demand that all controls yield the same dynamics in the absence of mutants: we set the wildtype mean and variance in this case to be identical under all controls. The mean is chosen to be *w*_{opt} (assuming that, without mutants, each control steers the population of wildtypes to its optimal value) and the variance is set by fixing the parameter *α*_{R} in the relaxed replication model to *α*_{R} = 10 (Eq (5)), as its value was originally estimated to lie in the range (5–17) [40]. This fixes the value for *c*_{1} in models (I)-(III) (S3 Table).

We now wish to optimise the parameters *δ* and *η* in models (III) and (IV). These parameters have direct impact on how mutants are sensed by the cell, and are important in determining how the dynamics change when mutant load increases. We thus want to investigate i) whether there exists an optimal amount of mutant sensing, and ii) how the cost of the dynamics resulting from these optimal parameters *δ*_{opt} and *η*_{opt} compares to that of our other two control strategies with *δ* fixed at either 0 or 1. To do this, we require both an optimization time-scale *T* and a set of initial conditions. We use *T* = ∞, corresponding to the steady state limit, and initial heteroplasmy values in the range *h*_{0} ∈ [0, 0.2]; we later consider finite values of *T*.

Through stochastic simulations, we find that i) a control lacking any mutant contribution shows an exponential increase in cost over time, and ii) effects of particular control strategies are more pronounced in low copy number cells (Table 1()) (Fig 3). The relaxed replication rate control and our linear feedback function behave very similarly when *δ* and *η* take their optimal values. Cost variances, as well as mutant and wildtype dynamics, are shown in S6 Fig. Model II shows an increase in mean cost over time while mean mutant and wildtype copy numbers remain constant (Fig 3 and S6 Fig)—this is due to increases in copy number variances as argued previously.

**A + B**: Here we show the mean cost (3 × 10^{4} repeats) for the following four controls: linear feedback controls λ(*w*, *m*) = *μ* + *c*_{1}(*w*_{opt} − (*w* + *δm*)) with I) *δ* = 0, II) *δ* = 1, and III) *δ* = *δ*_{opt}, and IV) the optimised ‘relaxed replication control’ [40, 41] (Eq (5)). Controls were initialised in steady state at *h*_{0} = 0.15. Both figures used the saturating output model; figures (A) and (B) correspond to low and high copy number cells, respectively. We used *ϵ*_{1} = 0.3; other control parameters used are specified in section 5.3 in S1 File. **C: A control based on sensing mitochondrial energy output is generally a good strategy**. This plot shows the optimal value of *δ* in our linear control as a function of *ϵ*_{1}, for the linear and saturating model and for both low (*h*_{0} = 0.1, solid line) and high (*h*_{0} = 0.8, dashed line) initial heteroplasmies. Here we used *T* = 100 and high copy number values for both models. Similar plots for *T* = 10^{4} are shown in S7 Fig, section 5.3 in S1 File. In the linear model *δ*_{opt} becomes negative for low *ϵ*_{1} values; as mutant copy number increases, a negative *δ* leads to an increase in wildtype to compensate for the deficient mutants.

We now investigate how the optimal value of mutant sensing for the linear control (*δ*_{opt}) depends on timescale *T*, initial heteroplasmy *h*_{0} and the ‘mutant pathology level’ described by *ϵ*_{1}. Here, we use the term ‘mutation pathology level’ to refer to a lower energy production rate of mutants due to a lower resource consumption rate, while ‘mutant sensing‘, as explained earlier, is used as a more general term. Intuitively, values of *ϵ*_{1} ≃ 1 have *δ*_{opt} ≈ 1: when wildtypes and mutants are equivalent, having a steady state with *w* + *m* = *w*_{opt} is desirable.

Values for *δ*_{opt} were found for the linear and saturating models, with low and high initial heteroplasmy values, for *T* = 100 days (Fig 3C). Having *δ* ≈ 1 means wildtypes and mutants are fed back similarly, whereas when *δ* ≪ 1 mutants are fed back less. For very deficient mutants (low *ϵ*_{1}), a low *δ*_{opt} ensures that wildtype copy number remains close to its optimal value to compensate for the mutants. Generally, as *ϵ*_{1} decreases, *δ*_{opt} decreases (Fig 3C). Similar results are found for longer timescales *T* (section 5.3 in S1 File).

If mitochondrial energy outputs are sensed, the quantity ‘*w* + *δm*’ represents the mitochondrial energy production rate (power production). In this case, a mutant with low *ϵ*_{1} produces less energy and is thus sensed less (low *δ*). The relation between *ϵ*_{1} and *δ* now obeys the optimal trend shown in Fig 3C. Therefore, control strategies based on the energy status of the cell can often outperform controls based on mtDNA copy number (which always have *δ* = 1) or sensing mtDNA mass (which would work well for deficient deletion mutants, but would be suboptimal for deficient point mutations) (Table 1()). Control based on copy number is preferred when a mutant is nearly as functional as a wildtype, in which case energy output and copy number are very much related. We have not used the expression for energy output in our cost function as a control strategy itself because claims based on the linear function *w* + *δm* are more general than one based on the details of our cost function.

#### Locally optimal control strategies map the control space of mtDNA populations.

Our cost function allows us to identify locally optimal controls: controls that, for each state (*w*, *m*), move the system in the direction of the largest decrease in cost. The resulting dynamics are shown in Fig 4. When heteroplasmy is high, the main priority is not always to decrease mutant copy number, but to increase wildtype copy number even if this means an increase in mutant load (region 2 in Fig 4A). Only after wildtype copy number has sufficiently increased should the focus be on decreasing *m*. At high copy numbers, the optimal dynamics are to decrease all mtDNA in an evenhanded manner (region 1) rather than decreasing *m* at a faster rate than *w*. For the saturating model, there is a divergence point in the space of local optimal strategies, reflecting the two local cost minima (high wildtype and high mutant) observed earlier (Fig 2). Hence, there are several regions of state space where even for pathological mutants, reduction of mutant mtDNA alone is not always the optimal control strategy (Table 1()). Finally, the less pathological the mutants become (e.g. Fig 4B), the more the locally optimal control starts to resemble a linear control. In the linear output model, the optimal control always shows linear behaviour (Fig 4C).

Streamplots in (*w*, *m*)-space show the dynamics resulting from a locally optimal control, for various parameters of *ϵ*_{1}. At each point, arrows show the direction corresponding to the largest decrease in cost. Regions are coloured according to the magnitude of the decrease in cost when moving in the optimal direction. Black arrows illustrate general trends in these regions. **A**: Region (1) shows that at high copy numbers, both mutant and wildtype mtDNAs should be decreased in an evenhanded manner; region (2) shows the possibility that the optimal control involves an increase in mutant copy number. **B**: A higher value for the parameter *ϵ*_{1} is used, meaning mutants are less pathological. **C**: the locally optimal control for the linear output model more closely resembles a linear control. In both (A) and (B) we see a divergence point (red asterisk) illustrating the fact that both high mutant and high wildtype states constitute local attractors of low cost (as in Fig 2).

#### A parameterised model of artificial mtDNA control for disease treatment.

In the previous section we identified locally optimal control strategies. Of course, these strategies may not be achievable by the cell (e.g. the cell may not be able to decouple biogenesis of wildtype and mutant mtDNA). However, such controls may still be possible through human intervention. This is why, in this section, we model recently developed genetic treatments to artificially control mtDNA populations. We then combine this treatment model with our linear feedback control λ(*w*, *m*) and our cost function.

Mitochondrially targeted zinc finger nucleases (mtZFNs) [35, 36, 60–62] are able to produce shifts in heteroplasmy by specifically cutting mutant mtDNA. To develop quantitative theory to understand and tune the effects of these interventions, we model nuclease transfection as inducing selective increases in mtDNA degradation, on the background of the linear cellular feedback control introduced earlier. Our transfection model contains three parameters describing strength (*I*_{0}), duration (*b*), and selectivity (*ξ*) of nuclease treatment (Methods). We assume that for every mutant that is cleaved by the endonucleases, *ξ* wildtypes are cleaved [36]. For example, when *ξ* = 1 there is no distinction between mutants and wildtypes, and when *ξ* = 0 there is no off-target cleavage.

We fit these treatment parameters, as well as our feedback control parameters *c*_{1} and *δ*, to recently obtained experimental data [36]. These data involve heteroplasmy and mtDNA copy number measurements during iterative treatments with mtZFNs of 80% heteroplasmic human osteosarcoma 143B cybrid cells. Four sequential cycles of transfection and recovery were performed, where each recovery period lasted 28 days [36]. We use the data provided in Ref. [36] as well as additional data from this reference which was not explicitly provided in the paper. For the current study, we have collected new data consisting of: i) measurements of total mtDNA copy number in pre-treatment cells (which are used as initial conditions in our inference model), and ii) measurements of mtZFN expression profiles (S9 Fig).

#### A Bayesian description of nuclease treatment.

We use Metropolis sampling to obtain posterior distributions of the parameters *I*_{0}, *b*, *ξ*, *c*_{1} and *δ* (Methods). Bayesian credible intervals for heteroplasmy and total copy number values during four consecutive rounds of treatment are shown in Fig 5A and 5B, illustrating the ability of this simple model to capture the dynamics resulting from nuclease activity. A periodicity of 28 days was imposed, representing the experimental protocol (Methods).

We used Metropolis sampling to fit our model parameters to recently obtained experimental data [36]. Solid black lines correspond to the maximum a posteriori (MAP) prediction and vertical dashed lines indicate transfection events (once every 28 days). **A**: Drawing from our posterior distributions (5 × 10^{4} samples), we show the mean and 50% and 95% credible intervals of our heteroplasmy dynamics predictions during four rounds of treatment and recovery. Deterministic simulations were used. Crosses indicate data points from Ref. citeGammage16Near. **B**: Similar to figure (A), but showing relative total mtDNA copy numbers. **C + D**: Heteroplasmy and copy number dynamics were measured during a single round of treatment and recovery in a setting in which the mtZFN concentration was reduced by incorporating hammerhead ribozymes in the mtZFN backbone [36]. The credible intervals shown were obtained by sampling from the posterior distributions of the parameters *I*_{0}, *b*, *c*_{1} and *δ* obtained using the data in figures (A) and (B), and using *ξ* ∼ 0.15 which represents the maximum likelihood estimate of *ξ* using this low mtZFN concentration data (Methods). Error bars in figure (D) show standard deviations of experimental measurements [36].

Our inference suggests the selectivity parameter *ξ* to lie in the range 0.6–0.8, indicating high levels of off-target cleavage (S10 Fig). This is not surprising given the large drop in total copy number (as low as ∼5% of initial values) combined with a modest shift in heteroplasmy (from 0.8 to ∼ 0.6) upon the first treatment. Supporting the high off-target cleavage, mtZFNs not targeted to any mtDNA sequence reduced copy numbers to 25% of their original values [36].

We now investigate whether our model can account for additional data (obtained in Ref. [36]) consisting of heteroplasmy and copy number measurements in a setting in which the concentration of mtZFNs is reduced by incorporating hammerhead ribozymes in the mtZFN backbone (for details, see Ref. [36]). A single round of treatment and recovery in this setting led to a large shift in heteroplasmy, from *h* ≈ 0.8 to *h* ≈ 0.2, and a drop in copy number similar to the previous setting in which mtZFN concentrations were higher (after 24 hours, mtDNA copy number dropped to ∼20% of its original value) [36]. These observations indicate that lower mtZFN concentrations lead to the treatment being more selective and, surprisingly, of similar strength. Because the additional data involves a different experimental setup inducing a large increase in selectivity, we adjust the parameter *ξ* to fit this additional data by finding its maximum likelihood estimate (Methods) but use posterior samples for all other parameters (obtained from inference based only on the data shown in Fig 5A and 5B). We find, consonant with an improved selectivity of this modified protocol, that *ξ* ≈ 0.15 in this low mtZFN concentration setting (Methods) and that our model can reproduce the heteroplasmy and copy number dynamics using our previously fitted parameters *I*_{0}, *b*, *c*_{1} and *δ* (Fig 5C and 5D).

Finally, we measured transient expression profiles of mtZFNs using the same transfection protocol as in Ref. [36] (S9 Fig). Posterior samples of the parameters *I*_{0} and *b* predict that mtZFN levels have dropped to very low levels 5 days post-transfection in the setting without hammerhead ribozymes, consistent with our obtained experimental data (S11 Fig). We thus show that our model is capable of capturing the dynamics of several data sets. Our mtZFN treatment model predicts that total copy number reaches a minimum at around 3 days (Fig 5B and 5D).

#### Knowledge of the heteroplasmy distribution of a tissue is important in determinining how effciently the tissue can be treated.

To explore the effect of the heteroplasmy distribution on treatment efficacy, we consider three initial *h* distributions with different variances but identical homogenate means. We treat these populations multiple times using the parameter fits obtained in the previous section. The resulting shifts in heteroplasmy distribution, including mean and threshold-crossing probability, are shown in Fig 6A and 6B. High heteroplasmy variances require many cells close to the two extremes *h* = 0 and *h* = 1, which are challenging to shift. A striking reduction in treatment efficacy is predicted as heteroplasmy variance increases while fixing its mean (Fig 6A and 6B). Threshold crossing probability (for example, *P*(*h* > 0.6)) also becomes harder to shift at higher variance. We conclude that tissues with a high mean heteroplasmy level will generally be harder to treat if the heteroplasmy variance is high, especially if this high mean level is caused by a small percentage of cells.

**A + B**: The effect of four simulated consecutive treatments on three different initial heteroplasmy distributions is shown; all initial distributions have identical means (〈*h*〉 = 0.8) but different variances (increasing from left to right). The higher the variance of the initial population, the harder to shift mean heteroplasmy values; mean values after each treatment as well as *P*(*h* > 0.6) are shown in figure (B). In these simulations we assumed that every cell gets transfected. **Gentle but sustained treatments induce larger heteroplasmy shifts than hard and brief treatments. C**: Both the linear and saturating model show a sharp drop in the optimal treatment strength *I*_{0,opt} as the mutants become more functional (i.e. as *ϵ*_{1} increases). **D**: Means and variances of mutant and wildtype copy numbers were simulated during a round of treatment and recovery, using: i) parameters fitted to the data shown in Fig 5A and 5B (blue), ii) a longer treatment duration (smaller *b*, green) and iii) a higher selectivity (smaller *ξ*, magenta). MtZFN levels first drop below 5% of their maximum values after ∼4.5 and ∼35 days for the short (blue, magenta) and long (green) treatments, respectively. The longer weaker treatment induces higher heteroplasmy shifts than the shorter stronger treatment. *ϵ*_{1} = 0.2 was used, the corresponding cost heatmap is shown. Error bars show standard deviations (based on 10^{4} stochastic simulations), further detailed are given in section 6.6 in S1 File. **E**: This figure also illustrates that gentle sustained treatments lead to larger heteroplasmy shifts. Examples of treatment trajectories are shown; after a single treatment, an initial heteroplasmy of 0.8 is mapped to 0.53 (short strong treatment) or 0.39 (long weak treatment). Parameters chosen in figures (A)–(E) are based on the inference performed earlier, their exact values are provided in section 6.6 in S1 File.

We can use our parameterised theory to find optimal treatment strengths *I*_{0,opt} for a given system. Fig 6C shows *I*_{0,opt} as a function of *ϵ*_{1}. Intuitively, the strongest treatment should be given to the least functional mutants, and when mutants are almost as functional as wildtypes it is preferable not to treat at all. The optimal treatment strength drops rather sharply as *ϵ*_{1} increases, and does so sooner for the saturating model. This last observation may be because at some point reducing heteroplasmy becomes more expensive as can be seen in Fig 2B. Optimal treatment strengths for longer treatments (higher *b*) show similar qualitative behaviour.

Fig 6D shows trajectories in (*w*, *m*) space throughout a single treatment and recovery phase. The three trajectories shown correspond to: i) a short and strong treatment, ii) a long and weak treatment, and iii) a short but more selective treatment. The value for *I*_{0} was chosen such that, for the specific treatment duration *b* used and given a fixed total simulation time, the shift in heteroplasmy was largest. For the short treatment a relatively large proportion of time is spent fluctuating around steady state values (dynamics which do not change mean *h*) due to a relatively quick recovery, whereas for longer treatments more time is spent in the treatment phase itself (dynamics which lower mean *h*). We thus find that in a given time frame, treating longer but weaker results in a lower final heteroplasmy value than treating short and strong (Table 1()). A weaker treatment also reduces the chance of a cell losing all its mtDNA molecules. Intuitively, a more selective treatment leads to larger heteroplasmy shifts.

Nuclease treatment and a subsequent recovery phase will have the net effect of mapping an initial heteroplasmy value *h*_{i} to a mean final heteroplasmy value, *h*_{f}. We simulated this mapping in the presence of cellular feedback control (Fig 6E), finding that heteroplasmy shifts are largest for intermediate heteroplasmies. The difference in treatment results for long compared to short treatments is also illustrated. Interestingly, for high *h* values, it is possible to end up with a **higher** heteroplasmy value after treatment, especially if *ξ* ≃ 1 (S8 Fig).

## Discussion

In this work, we have built a quantitative theory bridging stochastic optimal control, costs of mtDNA populations, and gene therapies. Our results contribute to a growing body of evidence [63–66] that the variance of mtDNA populations has important physiological and therapeutic implications independently of mean heteroplasmy, and underline that stochastic theory is required to understand this biologically and medically important quantity.

Key findings of our model (Table 1) include () the identification of tradeoffs in the control of one or the other mtDNA species; () the observation that increasing mtDNA variance can lead to increased energetic costs over time and ageing even when means and demands are preserved; () intermediate heteroplasmy states can be more expensive than states homoplasmic in either mutant or wildtype; () mutant sensing can be required to avoid an exponentially increasing cost; () sensing of cellular energetic status can be more effective than other targets like mitochondrial mass; () reduction of mutant mtDNA alone is not always the optimal control strategy; () high heteroplasmy variance challenges gene therapy treatments; and () weak, long gene therapy trajectories are more effective than short, intense ones.

Our findings hold qualitatively under the range of conditions we discuss above. The aim of our manuscript is not to make detailed quantitative predictions and conclusions based on complex models, nor do we intend to imply that our models are the only possible models one could construct. Rather, we aim to provide general biologically plausible models to gain qualitative insights and to comment on large-scale behaviours. To this end, our cost function, used to illustrate some of our results, is phenomenological and contains several parameters. Most of these are biologically interpretable, meaning their values can be obtained or estimated from the literature. The main elements in our cost function are quite general: terms involving supply, demand, and resource.

To test the qualitative shape of our cost function, one could sort cells based on mitochondrial copy number and heteroplasmy to obtain samples at different points in (*w*, *m*) space. Measurements of e.g. cell proliferation, ROS or apoptosis rates allow for the evaluation of an effective cost at each of these points. By measuring the relative consumption rates of NADH and succinate, as well as the amount of ATP produced per glucose consumed, in identical cells exposed to different energy demands, the saturating output model may be probed.

If the parameter *δ* is low, i.e. mutants are sensed less, mutant copy numbers at high heteroplasmies will be higher than wildtype copy numbers at low heteroplasmies. Experimentally, it has been observed that heteroplasmic cells can have total mtDNA copy number values that are 5-17-fold higher compared to cells homoplasmic in wildtype [67–70]. The cell has somehow allowed these mutants to expand, which may mean that they are less tightly controlled; controls based on total energy output or mtDNA mass (which can result in *δ* < 1) may lead to such behaviours. A control on mtDNA mass could explain why deletion mutants are often seen to expand [71, 72] and would also predict normal copy number levels in cells harbouring mtDNA point mutations. Recently, it was found that samples with mtDNA indels had very high mtDNA copy number levels, but single nucleotide variants did not [73].

We showed that heteroplasmy distributions in cell populations can provide important information about the possibility of successfully treating these cells with endonucleases. A tissue may be harder to treat if its high mean heteroplasmy level is caused by a small percentage of dysfunctional cells. Experimental values of mean homogenate heteroplasmy in heart tissue of patients with the 3243A>G mutation are roughly around 0.8 (though ranges can be large [74–77]) and muscle tissue often shows mosaic structures, with deficient patches of cells adjacent to healthy cells. These examples show that it may be that, at least in some cases, high mean levels are indeed caused by a relatively low percentage of cells, meaning that there are still challenges ahead for efficiently treating these tissues.

One of the features of our cost function is that resource limitations play an important role in shaping the cost landscape. There are indications that cellular levels of NAD (a coenzyme involved in oxidative phosphorylation) are limiting, and that a sufficient supply of NAD to mitochondria becomes critical [78–81]. An increase of intracellular NAD can lead to an increase in oxygen consumption and ATP production [81] indicating that resource limitation may, at least in some cases, be a genuine constraint. Adding various kinds of resources can significantly change mitochondrial basal respiration rate [82–84].

Like any other model, our models have a defined range of applicability. A key baseline assumption was using identical replication and degradation rates for mutants and wildtypes. Various possibilities of distinct rates have been offered in the literature, including faster mutant replication rates [22, 68, 85–88], lower mutant degradation rates [89], and higher mutant degradation rates [90, 91]. Including such differences, and other features such as *de novo* mutations, degradation control, and cell divisions [38, 64, 92, 93], constitute natural extensions to our theory.

## Methods

### Wildtype and mutant mtDNA evolution equations

Wildtype and mutant mtDNA copy numbers are considered to have birth rate λ(*w*, *m*) = *μ* + *c*_{1}(*w*_{opt} − (*w* + *δm*)) and death rate *μ*, leading to the following evolution equations:
(4)

The corresponding stochastic system, required to e.g. describe fixation, does not have an explicit solution due to nonlinearities. The deterministic steady state solution of Eq (4) is given by (*w*_{ss} + *δm*_{ss}) = *w*_{opt} and represents a straight line in (*w*, *m*)-space (S1A Fig), whose slope depends on the value of *δ*. Stochastic dynamics will fluctuate around the steady state line, causing heteroplasmy to change over time until fixation of either species occurs. This means that, over long times, a cell will reach either *h* = 0 or *h* = 1 (in the absence of mutations). When mutations do occur, a cell will always reach a state with *h* = 1 (though many different mutant species may be present).

### Relaxed replication model

The relaxed replication model assumes a constant death rate *μ* and a birth rate of the form
(5)
with *α*_{R} > 1 and *η* constants [40, 41]. We have renamed the parameters of the original model for convenience. Note that both *α*_{R} and *η* influence the mutant contribution to λ(*w*, *m*) (rather than the single parameter *δ* in our linear model).

### Expected cost per unit time

Let the cost per unit time of state (*w*, *m*) be denoted by *C*, and the cost corresponding to the steady state (*w*_{ss}, *m*_{ss}) by . Even if steady state copy numbers are constant over time (i.e. the mean values of *w* and *m* are always equal to *w*_{ss} and *m*_{ss}) the mean cost per unit time is generally not equal to . By performing a Taylor expansion, the mean cost per unit time can be written as follows:
(6)
where *E*[*C*](*t*) is the expected cost per unit time given that the trajectory starts in state (*w*_{ss}, *m*_{ss}), and all partial derivatives are evaluated at steady state. These findings imply the following: suppose all cells in a population of cells are initialised in a state with minimum cost (corresponding to some specific number of mutant and wildtype mtDNA molecules). At some later time, the mtDNA populations in the different cells will have drifted apart and even if mean copy numbers (averaged over all cells) of *w* and *m* are identical to their initial values, the increase in variance between cells means that the overall mean cost (averaged over all cells) is higher than it was initially.

### Cost function structure

We assume that the net energy supply per unit time in a state (*w*, *m*), called *S*(*w*, *m*), involves the following four terms: (i) the energy output per unit time (*s*_{i}) produced by the mitochondria; (ii) a maintenance cost per unit time (*ρ*_{1}) to maintain the mitochondria, as their presence imposes some energetic cost (e.g. mRNA and protein synthesis); (iii) a building cost (*ρ*_{2}) for the biogenesis of new mitochondria; and (iv) a degradation cost (*ρ*_{3}) to degrade mitochondria. We will assume that every mtDNA molecule is associated to a particular amount of mitochondrial volume which we refer to as a ‘mitochondrion’ (section 4 in S1 File).

At any time, mitochondria experience a certain energy demand and to meet this demand they need to have a certain resource consumption rate *r*_{i} (where *i* = *w*, *m* refers to wildtype or mutant). Here we use the term ‘resource’ as an amalgamation of the substrates used for the oxidation system. We need to specify the relationship between the power supply (*s*) and the rate of resources consumed (*r*_{i}) by mitochondria. We use two different models *s*(*r*_{i}) which are discussed further in section 3 in S1 File (7)
where *ϕ*, *β*, *k* and *s*_{max} are constants respectively describing the mitochondrial efficiency, a basal proton leak-like term, the saturation rate of the efficiency, and the maximum power supply (section 4 in S1 File).

We assume that pathological mutants can have a deficient electron transport chain (which may support a smaller flux leading to a lower resource consumption rate for mutants and therefore a lower ATP production rate) and a lower energy production efficiency, leading to the following mutant energy output: *ϵ*_{2}*s*(*ϵ*_{1}*r*_{w}). Here, *ϵ*_{1}, *ϵ*_{2} ∈ [0, 1] describe the mutant resource uptake rate and the mutant energy production efficiency relative to that of a wildtype, respectively. In the main text we set *ϵ*_{2} = 1; other values of *ϵ*_{2} are discussed in section 4.7 in S1 File.

The mitochondrial maintenance cost is denoted by *ρ*_{1} and corresponds to the energetic cost required to maintain the mitochondrion that contains the mtDNA. This energetic costs involves factors like the synthesis and degradation of mitochondrial proteins and enzymes. We assume the maintenance cost is the same for wildtype and mutant mitochondria (though for some mutations this is quite possibly not the case). The net energy supply per unit time, *S*(*w*, *m*), then follows as Eq 3.

To determine the value of *r*_{w} for a given state (*w*, *m*), we first check whether the demand *D* (which we assume is a constant) can be satisfied by supply *S*(*w*, *m*). If it can, we set Eq (3) equal to *D* and solve for *r*_{w}, i.e. we assume that if possible, the mitochondria will exactly satisfy demand. It may, however, not be possible to satisfy demand, which can be because of two reasons: i) there are not enough mitochondria present to produce enough energy, or ii) the resource supply rate, *R* (a constant), is not enough to meet demand. In the former case, we set *r*_{w} = *r*_{max} (a specified maximum resource consumption rate per mitochondrion): the mitochondria work as hard as possible to keep their energy output closest to demand. In the latter case, we assume that the total available resource supply is shared equally between the mitochondria: . Further details of the cost function are given in sections 3–5 in S1 File.

The parameters used in our cost function are summarised in S2 Table and motivated in section 4 in S1 File. Despite our model being simple, most parameters are biologically interpretable.

### Modelling control through mitochondrially targeted endonucleases

Experimentally, cells are transfected with two mtZFN monomers: one which binds selectively to mutant mtDNAs, and one that binds mutants and wildtypes with equal strength [62]. We simplify this picture by assuming an ‘effective’ mtZFN pool and use [*ZFN*] to denote its concentration. The increase in mtDNA degradation rate caused by the mtZFNs is then assumed to be proportional to [*ZFN*].

Nucleases are imported into the cell and then degrade over time, meaning that their concentration in the cell (and in the mitochondria) may be approximated by an immigration-death model:
(8)
where *I*(*t*) and *μ*_{Z} are the immigration and death rates of the effective mtZFN pool, respectively. In recent experiments [36], nucleases are expressed for short times meaning that the immigration rate will increase sharply at the start of the treatment after which it decreases over time: we chose to model *I*(*t*) as an exponentially decaying function, *I*(*t*) = *I*_{0}*e*^{−bt}, where *I*_{0} denotes the initial rate directly after the treatment is initiated and *b* is a constant describing the duration of the treatment. The mtZFN concentration now becomes
(9)
which is shown for various parameter values in S8A Fig. The data we use to fit our models concerns heteroplasmy and total copy number measurements over four rounds of treatment, each treatment consisting of mtZFN transfection followed by a 28-day recovery period. During this recovery period, total copy numbers recover their initial values due to cellular feedback control. The increase in mtDNA death rate due to the presence of the mtZFNs, *μ*_{ZFN}, is given by
(10)
where *i* = 0, 1, 2, 3 indicates the treatment round. This equation is simply stating that new mtZFNs are added every 28 days. Death rates for *m* and *w* are now assumed to be
(11)
where *μ* denotes the baseline degradation rate and *ξ* represents treatment selectivity (e.g. when *ξ* = 0 there is no off-target cleavage).

### Model fits using Metropolis sampling

To fit our nuclease model to recently obtained experimental data [36], we use Eq (4) with *μ* replaced by *μ*(*t*)_{w} or *μ*(*t*)_{m} and λ(*w*, *m*) given by Eq (1):
(12)

Total mtDNA copy numbers in pre-treatment 80% heteroplasmy cells were measured using quantitative PCR (section 6.4 in S1 File) and were found to be 889 ± 214 (S.E., *n* = 3). We therefore assume an initial total copy number of 900, meaning *w* and *m* were initialized at 0.2 ⋅ 900 = 180 and 0.8 ⋅ 900 = 720, respectively. These evolution equations incorporate cellular feedback control as well as the nuclease treatment which occurs in cycles of 28 days. The mtZFN degradation rate was assumed to be *μ*_{z} = ln(2) day^{−1}, corresponding to a half-life of 1 day. This is in accord with the experimental observation that almost no mtZFN was present 4 days post-transfection (with a half-life of 1 day, only 6% of initial copy numbers remain after 4 days).

MCMC inference was performed using the Python package Pymc3, a package designed for Bayesian statistical modelling and probabilistic machine learning [94]. A Gaussian error model was assumed, i.e. the observed heteroplasmy and total copy number data are given by
(13)
where and denote our predicted heteroplasmy and copy number values obtained by numerically solving Eq (12), and we allow for different noise variances for *h* and *T* (in general, different experimental errors are expected as different methods are used to measure *h* and *T*). A metropolis sampler is used for parameter estimation. Maximum a posteriori (MAP) values were found to be . Due to a degeneracy in our mtZFN dynamics model (section 6.5 in S1 File) the MAP values of *I*_{0} and *b* are not necessarily unique at large *b* (details in section 6.5 in S1 File).

We explore the ability of our model to account for additional data from Ref. [36] (Fig 5C and 5D) which was not included in our inference. Using the MAP values for parameters *I*_{0}, *b*, *c*_{1}, *δ*, and (based on the data shown in Fig 5A and 5B), the maximum likelihood estimate of *ξ* is obtained based on the additional data, using a Gaussian error model similar to Eq (13). This maximum likelihood value is *ξ* ≈ 0.15.

## Supporting information

### S1 Fig. A linear feedback control has straight steady state lines.

**A)** The deterministic steady state lines of the feedback control given in Eq 4, using our linear version of λ(*w*, *m*), are shown in (*w*, *m*) space for various values of *δ* (grey lines show particular examples of ranges of *δ*). Constant heteroplasmy lines form straight lines through the origin. **B, C, D) Equal variances for different feedback control mechanisms**. Three different controls (see legend), all of the form λ(*w* + *δm*) with *δ* = 0.5, show nearly identical wildtype, mutant and heteroplasmy variances. Other parameters used are *N*_{ss} = 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).

https://doi.org/10.1371/journal.pcbi.1007023.s001

(EPS)

### S2 Fig. 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) in S1 File. For the linear model (corresponding to the straight lines) the parameters *ϕ* and *β* are changed by 10%, for the saturating model we vary *s*_{max} and *k*. The magenta line indicates *r*_{max}. **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.

https://doi.org/10.1371/journal.pcbi.1007023.s002

(EPS)

### S3 Fig. Changing mutant efficiency (*ϵ*_{2}) does not lead to expensive intermediate heteroplasmies.

**A), B)** Similar to Fig 2 in the main text, these figures show the cost values in (*w*, *m*) space, but now as a function of *ϵ*_{2} (mutant efficiency) instead of *ϵ*_{1}. This time we show the cost in the entire space. The white lines show the region in which demand is satisfied for our default parameter values. Because mutants consume the same amount of resource as wildtypes (*ϵ*_{1} = 1), resource becomes limiting at relatively low values of *m* compared to when *ϵ*_{1} < 1. Note that intermediate heteroplasmies are not less efficient here.

https://doi.org/10.1371/journal.pcbi.1007023.s003

(EPS)

### S4 Fig. Intermediate h values require more resources to satisfy demand, but only if mutants consume less resources.

**A)** The resource consumption rates and energy production rates of wildtypes and mutants are shown for two states: (*w*_{1}, *m*_{1}, *h*_{1}) = (9000, 1000, 0.1) and (*w*_{2}, *m*_{2}, *h*_{2}) = (7000, 3000, 0.3). In both cases, the total energy output is equal to the demand. When heteroplasmy is higher (*h* = 0.3), the individual resource consumption rates are higher in order to maintain a constant total energy output. Overall, the state with *h* = 0.1 uses the least resources (Eq. (20) in S1 File). *ϵ*_{1} = 0.35 was used. **B)** This figure is similar to figure (D) but now the two states (*w*_{1}, *m*_{1}, *h*_{1}) = (3000, 7000, 0.7) and (*w*_{2}, *m*_{2}, *h*_{2}) = (1000, 9000, 0.9) are compared. The state with *h* = 0.9 uses the least resources (Eq. (21) in S1 File).

https://doi.org/10.1371/journal.pcbi.1007023.s004

(EPS)

### S5 Fig. The existence of intermediate heteroplasmy values is a robust feature of the saturating output model.

We show the value of *h*_{max}, the most expensive heteroplasmy value at constant copy number, as a function of total copy number and *ϵ*_{1} (describing mutant pathology). White regions correspond to *h*_{max} = 1. **A)** Using our default parameter values, an intermediate *h*_{max} exists for large enough mutant functionality *ϵ*_{1}. **B)** The parameter *s*_{max} 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 *s*_{max} is decreased by 50%. Again, change in *h*_{max} are small.

https://doi.org/10.1371/journal.pcbi.1007023.s005

(EPS)

### S6 Fig. Wildtype, mutant and cost dynamics for four different control strategies.

Dynamics are shown for the four controls I, II, III, and IV defined in Table 3 in S1 File. Again, we see that the effects of the control are more noticeable in low copy number cells. Parameters are set as given in Table 3 in S1 File. Values for *w*_{opt} 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.

https://doi.org/10.1371/journal.pcbi.1007023.s006

(EPS)

### S7 Fig. At long times and high heteroplasmies, energy sensing control becomes suboptimal.

The optimal value of *δ* in a linear feedback control is shown as a function of *ϵ*_{1}. Here we used *T* = 10^{4} days (optimization time) and low copy numbers for both the linear and saturating model. The solid and dashed lines correspond to trajectories starting at *h*_{0} = 0.1 and *h*_{0} = 0.8, respectively. The less resources the mutants consume (and the less output they therefore produce) the lower their optimal contribution to the control.

https://doi.org/10.1371/journal.pcbi.1007023.s007

(EPS)

### S8 Fig. Zinc finger nuclease concentrations for short and long treatments.

**A)** Here we show the concentration of mitochondrially targeted Zinc Fingers as modelled by Eq (9) in the main text. The parameter values for the short and strong treatment illustrated here (*I*_{0} = 36, *b* = 11) are similar to those found in fitting the model to the data. For the mtZFN degradation rate we used *μ*_{Z} = log(2) day^{−1} (corresponding to a mtZFN half-life of 1 day). **There exists a possibility of increasing heteroplasmy levels through treatment**. **B)** The probability of increasing heteroplasmy above its initial pre-treatment value *h*_{0}, after one round of treatment and recovery, is shown as a function of *h*_{0} and *ξ*. Cells are initialised with a total copy number of 500. The cross indicates the parameters used in figure (D). The parameter values for *I*_{0}, *b* and *c*_{1} are fixed at: (*I*_{0}, *b*, *c*_{1}) ≈ (39, 20, 3 × 10^{−4}); these values provide good fits to experimental data when assuming a total initial copy number of 500. We used *δ* = 1. **C)** Similar to figure (B), but now cells are initialised with a total copy number of 5000; in these large copy number cells stochastic fluctuations in copy number have less effect and the probabilities of exceeding initial heteroplasmy values are smaller compared to figure (B). **D)** An example of a distribution of post-treatment heteroplasmy values is shown using parameters *h*_{0} and *ξ* as indicated by the cross in figure (B). The orange line indicates the value of *h*_{0} (the heteroplasmy that was present before the treatment started).

https://doi.org/10.1371/journal.pcbi.1007023.s008

(EPS)

### S9 Fig. mtZFN expression profile during transient transfection of 143B cells.

**A)** Here we show a schematic of the experiments involving i) transient transfection of high-heteroplasmy cells with plasmids expressing mtZFN monomers and fluorescent marker proteins, ii) FACS-based selection of cells expressing both mtZFN monomers (NARPd(+) and COMPa(-)), and iii) phenotypic evaluation of treated cells. Technical details are provided in Ref. [35]. **B)** Western blots showing the mtZFN expression profile indicate that the mtZFNs are almost undetectable at 96 hours post-transfection, and completely undetectable at 120 hours. Details of the protocol are provided in Ref. [35].

https://doi.org/10.1371/journal.pcbi.1007023.s009

(EPS)

### S10 Fig. Posterior mtZFN treatment parameter distributions.

Here we show our posterior distributions obtained after running our MCMC algorithm (left) as well as the corresponding sample values (right). Prior distributions are provided in the text. The posterior of log_{10} *b* is cut off due to a degeneracy in our model (Eq. (22) in S1 File), which does not affect our model predictions.

https://doi.org/10.1371/journal.pcbi.1007023.s010

(EPS)

### S11 Fig. Predictions of mtZFN expression are broadly consonant with experimental data.

Drawing from our posterior distributions for *I*_{0} and *b* obtained through Metropolis sampling, we show 50% and 95% credible intervals of our predicted mtZFN expression profile (solid black line denotes the maximum a posteriori (MAP) estimate). Data points are obtained through quantification of the western blots shown in S9B Fig 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).

https://doi.org/10.1371/journal.pcbi.1007023.s011

(EPS)

### S1 Table. Analytical expressions for the means and variances according to the linear noise approximation.

Solutions are shown for wildtype, mutant, and heteroplasmy variances for various types of control. Dots indicate constant or exponentially decaying terms; full solutions are provided in Eqs. (11)–(13) in S1 File. Note that the initial rate of increase of heteroplasmy variance only depends on mtDNA copy number and turnover (see also [38]).

https://doi.org/10.1371/journal.pcbi.1007023.s012

(EPS)

### S2 Table. Parameters used in our mitochondrial cost function with their descriptions.

Parameter values are derived and motivated in Section 4 in S1 File. In this table, ‘high’ and ‘low’ refer to high and low copy numbers, respectively, and the abbreviations ‘sat.’ and ‘lin.’ are used to indicate our saturating and linear output models.

https://doi.org/10.1371/journal.pcbi.1007023.s013

(EPS)

### S3 Table. 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 [40] 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 *w*_{opt} are 1524 and 7616 (Table 2 in S1 File). We further use *μ* = 0.07 day^{−1}.

https://doi.org/10.1371/journal.pcbi.1007023.s014

(EPS)

## References

- 1.
Laar V. S. V. and Berman S. B. Mitochondrial dynamics in parkinson’s disease.
*Exp Neurol*218(2), 247–256 (2009). pmid:19332061 - 2.
Liesa M., Palacín M., and Zorzano A. Mitochondrial dynamics in mammalian health and disease.
*Physiol Rev*89(3), 799–845 (2009). pmid:19584314 - 3.
Grandemange S., Herzig S., and Martinou J.-C. Mitochondrial dynamics and cancer.
*Semin Cancer Biol*19(1), 50–56 (2009). pmid:19138741 - 4.
Zhu X., Perry G., Smith M. A., and Wang X. Abnormal mitochondrial dynamics in the pathogenesis of alzheimer’s disease.
*J Alzheimers Dis*33, S253–S262 (2013). pmid:22531428 - 5.
Chan D. C. Fusion and fission: Interlinked processes critical for mitochondrial health.
*Annu Rev Genet*46(1), 265–287 (2012). pmid:22934639 - 6.
Chen H. and Chan D. C. Mitochondrial dynamics—fusion, fission, movement, and mitophagy—in neurodegenerative diseases.
*Hum Mol Genet*18(R2), R169–R176 (2009). pmid:19808793 - 7.
Hayashi J.-L., Ohta S., Kikuchi A., Takemitsu M., Goto Y.-i., and Nonaka I. Introduction of disease-related mitochondrial DNA deletions into HeLa cells lacking mitochondrial DNA results in mitochondrial dysfunction.
*Proceedings of the National Academy of Sciences*88(23), 10614–10618 (1991). - 8.
Rossignol R., Faustin B., Rocher C., Malgat M., Mazat J.-P., and Letellier T. Mitochondrial threshold effects.
*Biochemical Journal*370(3), 751–762 (2003). pmid:12467494 - 9.
Attardi G., Yoneda M., and Chomyn A. Complementation and segregation behavior of disease-causing mitochondrial DNA mutations in cellular model systems.
*Biochimica et Biophysica Acta (BBA)-Molecular Basis of Disease*1271(1), 241–248 (1995). - 10.
Boulet L., Karpati G., and Shoubridge E. Distribution and threshold expression of the tRNA (Lys) mutation in skeletal muscle of patients with myoclonic epilepsy and ragged-red fibers (MERRF).
*American journal of human genetics*51(6), 1187 (1992). pmid:1334369 - 11.
Nakada K., Inoue K., Ono T., Isobe K., Ogura A., Goto Y.-i., Nonaka I., and Hayashi J.-I. Inter-mitochondrial complementation: mitochondria-specific system preventing mice from expression of disease phenotypes by mutant mtDNA.
*Nature medicine*7(8), 934–940 (2001). pmid:11479626 - 12.
Moraes C. T. and Schon E. A. [44] detection and analysis of mitochondrial DNA and RNA in muscle by in situ hybridization and single-fiber pcr.
*Methods in enzymology*264, 522–540 (1996). - 13.
Sciacco M., Bonilla E., Schon E. A., DiMauro S., and Moraes C. T. Distribution of wild-type and common deletion forms of mtDNA in normal and respiration-deficient muscle fibers from patients with mitochondrial myopathy.
*Human Molecular Genetics*3(1), 13–19 (1994). pmid:8162014 - 14.
Russell O. and Turnbull D. Mitochondrial DNA disease—molecular insights and potential routes to a cure.
*Experimental cell research*325(1), 38–43 (2014). pmid:24675282 - 15.
Durham S. E., Samuels D. C., Cree L. M., and Chinnery P. F. Normal levels of wild-type mitochondrial DNA maintain cytochrome c oxidase activity for two pathogenic mitochondrial DNA mutations but not for m. 3243a→ g.
*The American Journal of Human Genetics*81(1), 189–195 (2007). pmid:17564976 - 16.
Gitschlag B. L., Kirby C. S., Samuels D. C., Gangula R. D., Mallal S. A., and Patel M. R. Homeostatic responses regulate selfish mitochondrial genome dynamics in c. elegans.
*bioRxiv*, 050930 (2016). - 17.
Posakony J. W., England J. M., and Attardi G. Mitochondrial growth and division during the cell cycle in hela cells.
*The Journal of cell biology*74(2), 468–491 (1977). pmid:885911 - 18.
Jajoo R., Jung Y., Huh D., Viana M. P., Rafelski S. M., Springer M., and Paulsson J. Accurate concentration control of mitochondria and nucleoids.
*Science*351 (6269), 169–172 (2016). pmid:26744405 - 19.
Otten A. B., Theunissen T. E., Derhaag J. G., Lambrichs E. H., Boesten I. B., Winandy M., van Montfoort A. P., Tarbashevich K., Raz E., Gerards M., et al. Differences in strength and timing of the mtDNA bottleneck between zebrafish germline and non-germline cells.
*Cell Reports*16(3), 622–630 (2016). pmid:27373161 - 20.
Tang Y., Schon E. A., Wilichowski E., Vazquez-Memije M. E., Davidson E., and King M. P. Rearrangements of human mitochondrial DNA (mtDNA): new insights into the regulation of mtDNA copy number and gene expression.
*Molecular biology of the cell*11(4), 1471–1485 (2000). pmid:10749943 - 21.
Diaz F., Bayona-Bafaluy M. P., Rana M., Mora M., Hao H., and Moraes C. T. Human mitochondrial DNA with large deletions repopulates organelles faster than full-length genomes under relaxed copy number control.
*Nucleic acids research*30(21), 4626–4633 (2002). pmid:12409452 - 22.
Kowald A. and Kirkwood T. B. Transcription could be the key to the selection advantage of mitochondrial deletion mutants in aging.
*Proceedings of the National Academy of Sciences*111(8), 2972–2977 (2014). - 23.
Burgstaller J. P., Johnston I. G., and Poulton J. Mitochondrial dna disease and developmental implications for reproductive strategies.
*Molecular human reproduction*21(1), 11–22 (2015). pmid:25425607 - 24.
Røyrvik E., Burgstaller J. P., and Johnston I. G. mtDNA diversity in human populations highlights the merit of haplotype matching in gene therapies.
*Molecular Human Reproduction*22(11), 809–817 (2016). pmid:27609757 - 25.
Kim Y.-G., Cha J., and Chandrasegaran S. Hybrid restriction enzymes: zinc finger fusions to Fok i cleavage domain.
*Proceedings of the National Academy of Sciences*93(3), 1156–1160 (1996). - 26.
Urnov F. D., Miller J. C., Lee Y.-L., Beausejour C. M., Rock J. M., Augustus S., Jamieson A. C., Porteus M. H., Gregory P. D., and Holmes M. C. Highly efficient endogenous human gene correction using designed zinc-finger nucleases.
*Nature*435 (7042), 646–651 (2005). pmid:15806097 - 27.
Bibikova M., Golic M., Golic K. G., and Carroll D. Targeted chromosomal cleavage and mutagenesis in drosophila using zinc-finger nucleases.
*Genetics*161(3), 1169–1175 (2002). pmid:12136019 - 28.
Christian M., Cermak T., Doyle E. L., Schmidt C., Zhang F., Hummel A., Bogdanove A. J., and Voytas D. F. Targeting dna double-strand breaks with tal effector nucleases.
*Genetics*186(2), 757–761 (2010). pmid:20660643 - 29.
Li T., Huang S., Zhao X., Wright D. A., Carpenter S., Spalding M. H., Weeks D. P., and Yang B. Modularly assembled designer tal effector nucleases for targeted gene knockout and gene replacement in eukaryotes.
*Nucleic acids research*, gkr188 (2011). - 30.
Miller J. C., Tan S., Qiao G., Barlow K. A., Wang J., Xia D. F., Meng X., Paschon D. E., Leung E., Hinkley S. J., et al. A TALE nuclease architecture for efficient genome editing.
*Nature biotechnology*29(2), 143–148 (2011). pmid:21179091 - 31.
Joung J. K. and Sander J. D. TALENs: a widely applicable technology for targeted genome editing.
*Nature reviews Molecular cell biology*14(1), 49–55 (2013). pmid:23169466 - 32.
Bacman S. R., Williams S. L., Pinto M., Peralta S., and Moraes C. T. Specific elimination of mutant mitochondrial genomes in patient-derived cells by mitoTALENs.
*Nature Medicine*19(9), 1111–1113 (2013). pmid:23913125 - 33.
Reddy P., Ocampo A., Suzuki K., Luo J., Bacman S. R., Williams S. L., Sugawara A., Okamura D., Tsunekawa Y., Wu J., Lam D., Xiong X., Montserrat N., Esteban C. R., Liu G.-H., Sancho-Martinez I., Manau D., Civico S., Cardellach F., del Mar O’Callaghan M., Campistol J., Zhao H., Campistol J. M., Moraes C. T., and Belmonte J. C. I. Selective elimination of mitochondrial mutations in the germline by genome editing.
*Cell*161(3), 459–469 (2015). pmid:25910206 - 34.
Bacman S. R., Williams S. L., Pinto M., and Moraes C. T. The use of mitochondria-targeted endonucleases to manipulate mtDNA.
*Methods in enzymology*547, 373 (2014). pmid:25416366 - 35.
Gammage P. A., Van Haute L., and Minczuk M. Engineered mtZFNs for manipulation of human mitochondrial DNA heteroplasmy.
*Mitochondrial DNA: Methods and Protocols*, 145–162 (2016). - 36.
Gammage P. A., Gaude E., Van Haute L., Rebelo-Guiomar P., Jackson C. B., Rorbach J., Pekalski M. L., Robinson A. J., Charpentier M., Concordet J.-P., et al. Near-complete elimination of mutant mtDNA by iterative or dynamic dose-controlled treatment with mtZFNs.
*Nucleic Acids Research*, 7804–7816 (2016). pmid:27466392 - 37.
Hashimoto M., Bacman S. R., Peralta S., Falk M. J., Chomyn A., Chan D. C., Williams L. S., and Moraes C. T. Mitotalen: A general approach to reduce mutant mtDNA loads and restore oxidative phosphorylation function in mitochondrial diseases.
*Molecular Therapy*(2015). - 38.
Johnston I. G. and Jones N. S. Evolution of cell-to-cell variability in stochastic, controlled, heteroplasmic mtdna populations.
*The American Journal of Human Genetics*99(5), 1150–1162 (2016). pmid:27843124 - 39.
Gross N. J., Getz G. S., and Rabinowitz M. Apparent turnover of mitochondrial deoxyribonucleic acid and mitochondrial phospholipids in the tissues of the rat.
*Journal of Biological Chemistry*244(6), 1552–1562 (1969). pmid:5773057 - 40.
Chinnery P. F. and Samuels D. C. Relaxed replication of mtDNA: a model with implications for the expression of disease.
*The American Journal of Human Genetics*64(4), 1158–1165 (1999). pmid:10090901 - 41.
Capps G. J., Samuels D. C., and Chinnery P. F. A model of the nuclear control of mitochondrial DNA replication.
*Journal of theoretical biology*221(4), 565–583 (2003). pmid:12713941 - 42.
Elson J., Samuels D., Turnbull D., and Chinnery P. Random intracellular drift explains the clonal expansion of mitochondrial DNA mutations with age.
*The American Journal of Human Genetics*68(3), 802–806 (2001). pmid:11179029 - 43.
Kowald A. and Kirkwood T. B. Mitochondrial mutations and aging: random drift is insufficient to explain the accumulation of mitochondrial deletion mutants in short-lived animals.
*Aging cell*12(4), 728–731 (2013). pmid:23742009 - 44.
Lill R., Hoffmann B., Molik S., Pierik A. J., Rietzschel N., Stehling O., Uzarska M. A., Webert H., Wilbrecht C., and Mühlenhoff U. The role of mitochondria in cellular iron–sulfur protein biogenesis and iron metabolism.
*Biochimica et Biophysica Acta (BBA)-Molecular Cell Research*1823(9), 1491–1508 (2012). - 45.
Wang C. and Youle R. J. The role of mitochondria in apoptosis.
*Annual review of genetics*43, 95 (2009). pmid:19659442 - 46.
Rizzuto R., De Stefani D., Raffaello A., and Mammucari C. Mitochondria as sensors and regulators of calcium signalling.
*Nature reviews Molecular cell biology*13(9), 566–578 (2012). pmid:22850819 - 47.
Morán M., Moreno-Lastres D., Marín-Buera L., Arenas J., Martín M. A., and Ugalde C. Mitochondrial respiratory chain dysfunction: implications in neurodegeneration.
*Free Radical Biology and Medicine*53(3), 595–609 (2012). pmid:22595027 - 48.
Monternier P.-A., Marmillot V., Rouanet J.-L., and Roussel D. Mitochondrial phenotypic flexibility enhances energy savings during winter fast in king penguin chicks.
*Journal of Experimental Biology*217(15), 2691–2697 (2014). pmid:24803465 - 49.
Salin K., Roussel D., Rey B., and Voituron Y. David and Goliath: a mitochondrial coupling problem?
*Journal of Experimental Zoology Part A: Ecological Genetics and Physiology*317(5), 283–293 (2012). - 50.
Brand M. D., Harper M.-E., and Taylor H. C. Control of the effective P/O ratio of oxidative phosphorylation in liver mitochondria and hepatocytes.
*Biochemical Journal*291(3), 739–748 (1993). pmid:8489502 - 51.
Yadava N. and Nicholls D. G. Spare respiratory capacity rather than oxidative stress regulates glutamate excitotoxicity after partial respiratory inhibition of mitochondrial complex i with rotenone.
*The Journal of neuroscience*27(27), 7310–7317 (2007). pmid:17611283 - 52.
Hill B. G., Dranka B. P., Zou L., Chatham J. C., and Darley-Usmar V. M. Importance of the bioenergetic reserve capacity in response to cardiomyocyte stress induced by 4-hydroxynonenal.
*Biochemical Journal*424(1), 99–107 (2009). pmid:19740075 - 53.
Pfeiffer T., Schuster S., and Bonhoeffer S. Cooperation and competition in the evolution of atp-producing pathways.
*Science*292 (5516), 504–507 (2001). pmid:11283355 - 54.
Angulo-Brown F., Santillán M., and Calleja-Quevedo E. Thermodynamic optimality in some biochemical reactions.
*Il Nuovo Cimento D*17(1), 87–90 (1995). - 55.
Waddell T. G., Repovic P., Meléndez-Hevia E., Heinrich R., and Montero F. Optimization of glycolysis: new discussions.
*Biochemistry and Molecular Biology Education*27(1), 12–13 (1999). - 56.
van Dijken J. P., Weusthuis R. A., and Pronk J. T. Kinetics of growth and sugar consumption in yeasts.
*Antonie van leeuwenhoek*63(3-4), 343–352 (1993). pmid:8279829 - 57.
Sriskanthadevan S., Jeyaraju D. V., Chung T. E., Prabha S., Xu W., Skrtic M., Jhas B., Hurren R., Gronda M., Wang X., et al. Aml cells have low spare reserve capacity in their respiratory chain that renders them susceptible to oxidative metabolic stress.
*Blood*125(13), 2120–2130 (2015). pmid:25631767 - 58.
Pfleger J., He M., and Abdellatif M. Mitochondrial complex ii is a source of the reserve respiratory capacity that is regulated by metabolic sensors and promotes cell survival.
*Cell death & disease*6(7), e1835 (2015). - 59.
Sharpley M. S., Marciniak C., Eckel-Mahan K., McManus M., Crimi M., Waymire K., Lin C. S., Masubuchi S., Friend N., Koike M., et al. Heteroplasmy of mouse mtDNA is genetically unstable and results in altered behavior and cognition.
*Cell*151(2), 333–343 (2012). pmid:23063123 - 60.
Minczuk M., Papworth M. A., Kolasinska P., Murphy M. P., and Klug A. Sequence-specific modification of mitochondrial DNA using a chimeric zinc finger methylase.
*Proceedings of the National Academy of Sciences*103(52), 19689–19694 (2006). - 61.
Minczuk M., Papworth M. A., Miller J. C., Murphy M. P., and Klug A. Development of a single-chain, quasi-dimeric zinc-finger nuclease for the selective degradation of mutated human mitochondrial DNA.
*Nucleic acids research*36(12), 3926–3938 (2008). pmid:18511461 - 62.
Gammage P. A., Rorbach J., Vincent A. I., Rebar E. J., and Minczuk M. Mitochondrially targeted ZFNs for selective degradation of pathogenic mitochondrial genomes bearing large-scale deletions or point mutations.
*EMBO molecular medicine*6(4), 458–466 (2014). pmid:24567072 - 63.
Cree L. M., Samuels D. C., de Sousa Lopes S. C., Rajasimha H. K., Wonnapinij P., Mann J. R., Dahl H.-H. M., and Chinnery P. F. A reduction of mitochondrial dna molecules during embryogenesis explains the rapid segregation of genotypes.
*Nature genetics*40(2), 249–254 (2008). pmid:18223651 - 64.
Johnston I. G., Burgstaller J. P., Havlicek V., Kolbe T., Rülicke T., Brem G., Poulton J., and Jones N. S. Stochastic modelling, bayesian inference, and new in vivo measurements elucidate the debated mtDNA bottleneck mechanism.
*eLife*4, e07464 (2015). pmid:26035426 - 65.
Freyer C., Cree L. M., Mourier A., Stewart J. B., Koolmeister C., Milenkovic D., Wai T., Floros V. I., Hagström E., Chatzidaki E. E., et al. Variation in germline mtdna heteroplasmy is determined prenatally but modified during subsequent transmission.
*Nature genetics*44(11), 1282–1285 (2012). pmid:23042113 - 66.
Stewart J. B. and Chinnery P. F. The dynamics of mitochondrial dna heteroplasmy: implications for human health and disease.
*Nature Reviews Genetics*16(9), 530–542 (2015). pmid:26281784 - 67.
Kaufmann P., Shanske S., Hirano M., DiMauro S., King M. P., Koga Y., and Schon E. A. Mitochondrial DNA and RNA processing in melas.
*Annals of neurology*40(2), 172–180 (1996). pmid:8773598 - 68.
Shoubridge E. A., Karpati G., and Hastings K. E. Deletion mutants are functionally dominant over wild-type mitochondrial genomes in skeletal muscle fiber segments in mitochondrial disease.
*Cell*62(1), 43–49 (1990). pmid:2163769 - 69.
Moraes C. T., Ricci E., Petruzzella V., Shanske S., DiMauro S., Schon E. A., and Bonilla E. Molecular analysis of the muscle pathology associated with mitochondrial DNA deletions.
*Nature genetics*1(5), 359–367 (1992). pmid:1284549 - 70.
Tokunaga M., Mita S., Murakami T., Kumamoto T., Uchino M., Nonaka I., and Ando M. Single muscle fiber analysis of mitochondrial myopathy, encephalopathy, lactic acidosis, and stroke-like episodes (melas).
*Annals of neurology*35(4), 413–419 (1994). pmid:8154867 - 71.
Nicholas A., Kraytsberg Y., Guo X., and Khrapko K. On the timing and the extent of clonal expansion of mtDNA deletions: evidence from single-molecule pcr.
*Experimental neurology*218(2), 316–319 (2009). pmid:19426731 - 72.
Schwarze S. R., Lee C. M., Chung S. S., Roecker E. B., Weindruch R., and Aiken J. M. High levels of mitochondrial DNA deletions in skeletal muscle of old rhesus monkeys.
*Mechanisms of ageing and development*83(2), 91–101 (1995). pmid:8569289 - 73.
Reznik E., Miller M. L., Şenbabaoğlu Y., Riaz N., Sarungbam J., Tickoo S. K., Al-Ahmadie H. A., Lee W., Seshan V. E., Hakimi A. A., et al. Mitochondrial DNA copy number variation across human cancers.
*eLife*5, e10769 (2016). pmid:26901439 - 74.
Majamaa-Voltti K., Peuhkurinen K., Kortelainen M.-L., Hassinen I. E., and Majamaa K. Cardiac abnormalities in patients with mitochondrial DNA mutation 3243A>G.
*BMC cardiovascular disorders*2(1), 12 (2002). pmid:12150714 - 75.
Ciafaloni E., Ricci E., Servidei S., Shanske S., Silvestri G., Manfredi G., Schon E., and DiMauro S. Widespread tissue distribution of a trnaleu (uur) mutation in the mitochondrial DNA of a patient with melas syndrome.
*Neurology*41(10), 1663–1663 (1991). pmid:1922812 - 76.
Anan R., Nakagawa M., Miyata M., Higuchi I., Nakao S., Suehara M., Osame M., and Tanaka H. Cardiac involvement in mitochondrial diseases a study on 17 patients with documented mitochondrial DNA defects.
*Circulation*91(4), 955–961 (1995). pmid:7850981 - 77.
Jeppesen T. D., Schwartz M., Olsen D. B., and Vissing J. Oxidative capacity correlates with muscle mutation load in mitochondrial myopathy.
*Annals of Neurology*54(1), 86–92 (2003). pmid:12838523 - 78.
Stein L. R. and Imai S.-i. The dynamic regulation of nad metabolism in mitochondria.
*Trends in Endocrinology & Metabolism*23(9), 420–428 (2012). - 79.
Alano C. C., Ying W., and Swanson R. A. Poly (adp-ribose) polymerase-1-mediated cell death in astrocytes requires NAD+ depletion and mitochondrial permeability transition.
*Journal of Biological Chemistry*279(18), 18895–18902 (2004). pmid:14960594 - 80.
Bai P., Cantó C., Oudart H., Brunyánszki A., Cen Y., Thomas C., Yamamoto H., Huber A., Kiss B., Houtkooper R. H., et al. PARP-1 inhibition increases mitochondrial metabolism through SIRT1 activation.
*Cell Metabolism*13(4), 461–468 (2011). pmid:21459330 - 81.
Pittelli M., Felici R., Pitozzi V., Giovannelli L., Bigagli E., Cialdai F., Romano G., Moroni F., and Chiarugi A. Pharmacological effects of exogenous nad on mitochondrial bioenergetics, DNA repair, and apoptosis.
*Molecular pharmacology*80(6), 1136–1146 (2011). pmid:21917911 - 82.
Korzeniewski B., Harper M.-E., and Brand M. D. Proportional activation coefficients during stimulation of oxidative phosphorylation by lactate and pyruvate or by vasopressin.
*Biochimica et Biophysica Acta (BBA)-Bioenergetics*1229(3), 315–322 (1995). - 83.
Nobes C., Hay W., and Brand M. The mechanism of stimulation of respiration by fatty acids in isolated hepatocytes.
*Journal of Biological Chemistry*265(22), 12910–12915 (1990). pmid:2376580 - 84.
Affourtit C. and Brand M. D. Uncoupling protein-2 contributes significantly to high mitochondrial proton leak in INS-1E insulinoma cells and attenuates glucose-stimulated insulin secretion.
*Biochemical Journal*409(1), 199–204 (2008). pmid:17824844 - 85.
Michikawa Y., Mazzucchelli F., Bresolin N., Scarlato G., and Attardi G. Aging-dependent large accumulation of point mutations in the human mtDNA control region for replication.
*Science*286 (5440), 774–779 (1999). pmid:10531063 - 86.
Wallace D. C.
Mitochondrial genetics: a paradigm for aging and degenerative diseases?
*Science*256 (5057), 628 (1992). - 87.
Kowald A. and Klipp E. Mathematical models of mitochondrial aging and dynamics.
*Prog Mol Biol Transl Sci*127, 63–92 (2014). pmid:25149214 - 88.
Yoneda M., Chomyn A., Martinuzzi A., Hurko O., and Attardi G. Marked replicative advantage of human mtDNA carrying a point mutation that causes the melas encephalomyopathy.
*Proceedings of the National Academy of Sciences*89(23), 11164–11168 (1992). - 89.
Kowald A. and Kirkwood T. B. Accumulation of defective mitochondria through delayed degradation of damaged organelles and its possible role in the ageing of post-mitotic and dividing cells.
*Journal of theoretical biology*202(2), 145–160 (2000). pmid:10640434 - 90.
Twig G., Elorza A., Molina A. J., Mohamed H., Wikstrom J. D., Walzer G., Stiles L., Haigh S. E., Katz S., Las G., et al. Fission and selective fusion govern mitochondrial segregation and elimination by autophagy.
*The EMBO journal*27(2), 433–446 (2008). pmid:18200046 - 91.
Busch K. B., Kowald A., and Spelbrink J. N. Quality matters: how does mitochondrial network dynamics and quality control impact on mtDNA integrity?
*Philosophical Transactions of the Royal Society B: Biological Sciences*369 (1646), 20130442 (2014). - 92.
Poovathingal S. K., Gruber J., Halliwell B., and Gunawan R. Stochastic drift in mitochondrial DNA point mutations: a novel perspective ex silico.
*PLoS computational biology*5(11), e1000572 (2009). pmid:19936024 - 93.
Johnston I. G. and Jones N. S. Closed-form stochastic solutions for non-equilibrium dynamics and inheritance of cellular components over many cell divisions. In
*Proc. R. Soc. A*, volume 471, 20150050. The Royal Society, (2015). - 94.
Salvatier J., Wiecki T. V., and Fonnesbeck C. Probabilistic programming in python using pymc3.
*PeerJ Computer Science*2, e55 (2016).