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.
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.
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: λ (w, m) = μ + c 1 (w o p t - (w + δ m) ) (1) where c1 > 0, wopt > 0 and δ are constants, with wopt denoting the steady state value towards which the effective population, here defined as w + δm, is controlled. The magnitude of c1 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 = wopt), 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 wopt 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 wopt. 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. 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.
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 (I) the identification of tradeoffs in the control of one or the other mtDNA species; (II) the observation that increasing mtDNA variance can lead to increased energetic costs over time and ageing even when means and demands are preserved; (III) intermediate heteroplasmy states can be more expensive than states homoplasmic in either mutant or wildtype; (IV) mutant sensing can be required to avoid an exponentially increasing cost; (V) sensing of cellular energetic status can be more effective than other targets like mitochondrial mass; (VI) reduction of mutant mtDNA alone is not always the optimal control strategy; (VII) high heteroplasmy variance challenges gene therapy treatments; and (VIII) 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.
Wildtype and mutant mtDNA copy numbers are considered to have birth rate λ (w, m) = μ + c1 (wopt − (w + δm) ) and death rate μ, leading to the following evolution equations: d w d t = w (λ (w, m) - μ) d m d t = m (λ (w, m) - μ) (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 (wss + δmss) = wopt 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). The relaxed replication model assumes a constant death rate μ and a birth rate of the form λ (w, m) = μ w + m (α R [ w o p t - (w + η m) ] + w + η m) (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). Let the cost per unit time of state (w, m) be denoted by C, and the cost corresponding to the steady state (wss, mss) by C ¯. Even if steady state copy numbers are constant over time (i. e. the mean values of w and m are always equal to wss and mss) the mean cost per unit time is generally not equal to C ¯. By performing a Taylor expansion, the mean cost per unit time can be written as follows: E [ C ] (t) ≈ C ¯+ 1 2 (var (w (t) ) ∂ 2 C ∂ w 2 + var (m (t) ) ∂ 2 C ∂ m 2 + 2 cov (w (t), m (t) ) ∂ 2 C ∂ w ∂ m) (6) where E[C] (t) is the expected cost per unit time given that the trajectory starts in state (wss, mss), 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. 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 (si) 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 ri (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 (ri) by mitochondria. We use two different models s (ri) which are discussed further in section 3 in S1 File s (r w) = ϕ (r w - β) s (r w) = 2 s m a x 1 + e - k r w - 1. 1 s m a x (7) where ϕ, β, k and smax 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: ϵ2s (ϵ1rw). 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 rw 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 rw, 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 rw = rmax (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: r w = R w + ϵ 1 m. 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. 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: d [ Z F N ] (t) d t = I (t) - μ z [ Z F N ] (t) (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) = I0e−bt, where I0 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 [ Z F N ] (t) = I 0 μ z - b (e - b t - e - μ z t) (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 μ Z F N (28 · i < t < 28 · (i + 1) ) = μ + ∑ j = 0 i [ Z F N ] (t - 28 · j) (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 μ (t) w = μ + ξ · μ Z F N (t) μ (t) m = μ + μ Z F N (t) (11) where μ denotes the baseline degradation rate and ξ represents treatment selectivity (e. g. when ξ = 0 there is no off-target cleavage). 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): d w d t = w [ c 1 (w o p t - (w + δ m) ) - ξ · μ Z F N (t) ] d m d t = m [ c 1 (w o p t - (w + δ m) ) - μ Z F N (t) ] (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 y i (h) and total copy number y i (T) data are given by y i (h) = y ^ i (h) + N (0, σ h 2) y i (T) = y ^ i (T) + N (0, σ T 2) (13) where y ^ i (h) and y ^ i (T) 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 (I 0, b, c 1, ξ, δ, σ h 2, σ T 2) M A P ≈ (122. 82,46. 68,1. 90 × 10 - 4,0. 72,1. 26,0. 061,0. 10). Due to a degeneracy in our mtZFN dynamics model (section 6. 5 in S1 File) the MAP values of I0 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 I0, b, c1, δ, σ h 2 and σ T 2 (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.
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.
genome engineering
control theory
medicine and health sciences
gene therapy
mitochondrial dna
engineering and technology
nucleases
enzymes
synthetic biology
dna-binding proteins
enzymology
plant energy production
synthetic bioengineering
control engineering
plant science
systems science
mathematics
forms of dna
heteroplasmy
dna
bioenergetics
mitochondria
molecular biology techniques
cellular structures and organelles
synthetic genomics
research and analysis methods
bioengineering
synthetic genome editing
computer and information sciences
proteins
molecular biology
clinical genetics
biochemistry
plant biochemistry
hydrolases
zinc finger nucleases
cell biology
nucleic acids
heredity
genetics
biology and life sciences
physical sciences
energy-producing organelles
Energetic costs of cellular and therapeutic control of stochastic mitochondrial DNA populations