Closed-form stochastic solutions for non-equilibrium dynamics and inheritance of cellular components over many cell divisions

Stochastic dynamics govern many important processes in cellular biology, and an underlying theoretical approach describing these dynamics is desirable to address a wealth of questions in biology and medicine. Mathematical tools exist for treating several important examples of these stochastic processes, most notably gene expression and random partitioning at single-cell divisions or after a steady state has been reached. Comparatively little work exists exploring different and specific ways that repeated cell divisions can lead to stochastic inheritance of unequilibrated cellular populations. Here we introduce a mathematical formalism to describe cellular agents that are subject to random creation, replication and/or degradation, and are inherited according to a range of random dynamics at cell divisions. We obtain closed-form generating functions describing systems at any time after any number of cell divisions for binomial partitioning and divisions provoking a deterministic or random, subtractive or additive change in copy number, and show that these solutions agree exactly with stochastic simulation. We apply this general formalism to several example problems involving the dynamics of mitochondrial DNA during development and organismal lifetimes.

In this Appendix we present some of the lengthier mathematical results used in the main text. We include a Mathematica notebook as part of this article, containing the following derivations.

Generating function for the BID process
The derivation of the generating function for BID dynamics in the absence of cell divisions is well known: we include it here for completeness. We write the PDE describing the generating function G(z, t) in Laplace form: We proceed by using the method of characteristics, writing down ODEs describing how the parameters of G, and G itself, changes along a characteristic curve, with progress along such a curve parameterised by s. The corresponding ODEs are Eqn. 4 lets us immediately set t = s, omitting a constant of integration as the absolute value of progress along a characteristic curve is unimportant. Using t = s throughout, Eqn. 5 is solved by where c 1 is a constant of integration, the explicit form of which will be useful later. Rearranging this into an expression for c 1 gives Finally, Eqn. 6 with Eqn. 7 gives us G = c 2 e −αt e λt+c1ν − λe λc1+νt α/λ .
c 2 is a function of c 1 because the quantity c 1 , the constant of integration acquired when integrating z with respect to s, is independent of s, and hence forms an independent parameter when integrating G with respect to s. We require that G(t = 0) = z m0 , so we choose where the first term cancels the final term in Eqn. 9 when t = 0, and the final term can be seen to extract a factor z m0 from Eqn. 8 for c 1 when t = 0. We then have G(z, t) = c 2 (c 1 (z, t))e −αt e λt+c1(z,t)ν − λe λc1(z,t)+νt α/λ which, after inserting Eqn. 8 and some algebra, gives Recurrence relations arising from induction over cell divisions This section focusses on the solution of recurrence relations of the form In the Main Text, both h i and z i follow relationships of this kind; we use the symbol ζ i here to emphasise that the same solution strategy applies in both cases, and describe specific solutions below. This system is solved, after [1], by defining α ≡ a+d c , β ≡ D c 2 , D ≡ ad − bc, y i = ζ i + d c and implicitly defining w i through y i = wi wi+1 . These changes of variables allow us to find an expression for w i , which can then be substituted back through the above chain to find ζ i . We have which is solved by considering solutions to the characteristic equation βk 2 − αk + 1 = 0, which are straightforwardly k 1,2 = 1 2β (α ± α 2 − 4β). Then where C i are constants to be determined from boundary conditions. If the boundary condition takes the form ζ n = pz+q rz+s , as is the case throughout the situations we consider, we obtain Thus, given knowledge of a, b, c, d from the recurrence relation and p, q, r, s from the initial condition, we can obtain k 1 , k 2 through α and β and hence use Eqns. 22 and 20 to obtain a solution to the recurrence relation. Below, we use this approach to obtain solutions for the systems of interest in the main text.

Products of prefactors over the inductive process
We are concerned with an expression for the product n i=1 φ(z i , τ ) from the Main Text. We first consider n i=1 ξ(z i , τ ) which occurs as a factor in this expression in every partitioning regime. We recall that ξ(z, t) has the form If we write z i in the form a general form separating variables raised to the power i in z i (here represented by ρ), from their coefficientsÃ j ,B j , the form of Eqn. 45 gives us: where The product of an expression of this form can be written Here we make use of the q-Pochhammer symbol (a; q) n , defined by The terms within the brackets on the RHS of Eqn. 48 then become, by setting a = −A j /B j and q = ρ A /ρ B , and so yielding a simple form for the product of interest. For the binomial case, we have from Eqn. 38: Here, the two quantities raised to the power of i in z i are l and 2, so we set ρ A ≡ l, ρ b ≡ 2. Then by comparing For the subtractive case, we have from Eqn. 44: Now the quantities raised to the power of i are l and 1, so we set ρ A = l, ρ B = 1, and findÃ

Other products
We can also use this result to compute the product of exponentiated prefactors involved in the subtractive inheritance regimes. We will first consider the product n i=1 g(z i , τ ) −η , which plays a role in the deterministic subtractive inheritance we consider later. We recall the definitions of the recurrence relations in the Main Text For both subtractive inheritance cases, θ(g(z, t)) = g(z, t), so it can straightforwardly be seen that Thus, the product As i and j are dummy variables, we can then identify the required solution as n i=1 h −η i . We have from Eqn. 32 that where we have rewritten the final line to avoid a diverging factor of (z − 1) −1 appearing in the Pochhammer symbol, using We then see that h −η i is of the form Eqn. 47, so we can use the result therein, with Next, we wish to compute an expression for , of use in the random subtractive regime. We again exploit the relation between g(z i , τ ) and h j in Eqn. 59 to show that the kernel of the desired product is equivalent to 1 2 + 1 2hi . Again, we use h i from Eqn. 32; after some algebra this expression reduces to whereupon we can use Eqn. 52 with

Full forms of generating functions
To recap, we use α, λ, ν to respectively represent the rates of immigration, birth, and death in our model; m 0 for initial copy number; τ for cell cycle length, n for the number of divisions that have occurred, and t for the elapsed time since the most recent cell division. We employ simplifying symbols l ≡ e (λ−ν)τ ; l ≡ e (λ−ν)t and The general form of the generating functions was shown in the Main Text to be where z i and h i are the solutions to recursion relations defined in the Main Text. The term (iii) is the same in all calculations and is, from Eqn. 12, The generating function in the case of binomial partitioning at cell divisions involves (i) Eqn. 52 applied to the appropriate terms described in Eqn. 53; (ii) unity; (iii) Eqn. 66; and (iv) h 0 from Eqn. 27, giving overall The generating function in the case of random subtractive inheritance involves (i) Eqn. 52 applied to the appropriate terms described in Eqn. 54; (ii) Eqn. 52 applied to Eqn. 64 as described; (iii) Eqn. 66; and (iv) h 0 from Eqn. 32, and is

Birth-death dynamics for the mtDNA bottleneck
For birth-death dynamics, α = 0, so for binomial partitioning the generating function takes the form of the final term in Eqn. 67.

(72)
Relaxed replication of mtDNA dynamics The generating function for a single cell cycle (involving immigration and death dynamics) can straightforwardly be found and so For binomial partitioning, we have φ(z, t) = 1 and θ(g(z, t)) = (1/2 + g(z, t)/2). The solutions of the appropriate recurrence relations are then: So the overall generating function is which after some simplification can be written as In particular, the mean copy number is ae az+b (cz + d) m0 + e az+b m 0 c(cz + d) m0−1 z=1 , giving E(m, t) = 2 −n e −β(t+nτ ) m 0 + m opt 1 − e −βt + 2 −n e −βt (e βτ − 1)(2 n − e −βnτ ) 2e βτ − 1 (84) and setting t = τ (at the end of a cell cycle), we obtain The n → ∞ limit of this expression is In this n → ∞ limit the generating function reduces to so that Using Eqn. 79 and Leibniz's rule, the general probability distribution function is given by where U (a, b, z) is the confluent hypergeometric function. For subtractive partitioning, we have φ(z, t) = (1/2 + 1/(2g(z, t))) 2η and θ(g(z, t)) = g(z, t). The solutions of the appropriate recurrence relations are then: So the overall generating function is In particular, E(m, t) = m opt + e −β(t+nτ ) (m 0 − m opt + η(1 + (0; e βτ ) n+1 )).
It follows straightforwardly from the definition of the q-Pochhammer symbol that Using these results and setting t = τ (at the end of a cell cycle), we obtain after some manipulation and the only term retained in the n → ∞ limit is
In the absence of cell divisions, the birth-death generating function is simply Eqn. 12 with α = 0. Setting z = 0 gives the general extinction probability Copy number balance can be enforced in the absence of cell divisions by taking the λ = ν limit, from which follows the generating function from which straightforwardly follows the extinction probability P BD ,H (m = 0) = νt 1 + νt m0 (120)

Other partitioning regimes
We have derived results for the case where a binomially-distributed random number of agents is lost at each cell division. We now consider the case where this number is a fixed constant. We will denote this constant loss number by η. In this case, and so φ(z, t) = g(z, t) −η ; (123) θ(g(z, t)) = g(z, t).
As θ, ξ and g take the same form as for the random loss case, the solutions for z i and h i are the same as before. The difference (due to the different form of φ) is in the second product in the general generating function form, which is now n i=1 g(z i , τ ) −η . By Eqn. 61 we have that this factor takes the form of Eqn. 52, with The generating function in the case of deterministic subtractive inheritance involves (i) Eqn. 52 applied to the appropriate terms described in Eqn. 54; (ii) Eqn. 52 applied to Eqn. 61 as described; (iii) Eqn. 66; and (iv) h 0 from Eqn. 32, and is Now we briefly explore two other inheritance regimes of potential biological applicability. In these cases we have not been able to obtain closed-form solutions for an arbitrarily large number of cell divisions: however, the appropriate recursion relations may be followed for as many divisions as required in order to obtain a closed-form solution for the generating function.
Next we consider the binomial inheritance of clusters of agents. We will assume that these clusters are of fixed size n c . In this case, we consider the new variables C b = m i,b /n c , C a = m i,a /n c (denoting the number of clusters before and after a cell division), and write g mi,a P δ (m i,a |m i,b )P i−1 (m i,b , τ |m 0 ) The resultant generating function analysis yields a very similar outcome to those previously considered, with the altered recurrence relation z i = 1 2 + g(z i−1 , τ ) nc 2 1 nc ; z 1 = 1 2 + g(z, t) nc 2 1 nc (135) h i = g 0 1 2 + h nc i+1 2 1 nc , τ ; h n = g 0 (z, t).
We have been unable to reduce the recurrence relations Eqns. 130-131 or Eqns. 135-136 to a closed-form solution for birth-death dynamics, but the corresponding problems may be solved for an arbitrary number of cell divisions by writing out the recurrence explicitly, thereby obtaining the generating function for a given number of cell divisions. The figure in the main text uses this approach.