∎

11institutetext: Timothy C. Stutz 22institutetext: Department of Computational
Medicine
University of California, Los Angeles, CA
22email<EMAIL_ADDRESS>33institutetext: Janet S. Sinsheimer
44institutetext: Departments of Biostatistics, Computational Medicine, Human
Genetics
University of California, Los Angeles, CA
44email<EMAIL_ADDRESS>55institutetext: Mary E. Sehl 66institutetext:
Department of Computational Medicine and Division of Hematology-Oncology,
Department of Medicine
David Geffen School of Medicine, University of California, Los Angeles, CA
66email<EMAIL_ADDRESS>77institutetext: Jason Xu 88institutetext:
Department of Statistical Science
Duke University, Durham, NC
88email<EMAIL_ADDRESS>
corresponding author

# Computational tools for assessing gene therapy under branching process
models of mutation

Timothy C. Stutz  Janet S. Sinsheimer  Mary Sehl  Jason Xu

(Received: date / Accepted: date)

###### Abstract

Multitype branching processes are ideal for studying the population dynamics
of stem cell populations undergoing mutation accumulation over the years
following transplant. In such stochastic models, several quantities are of
clinical interest as insertional mutagenesis carries the potential threat of
leukemogenesis following gene therapy with autologous stem cell
transplantation. In this paper, we develop a three-type branching process
model describing accumulations of mutations in a population of stem cells
distinguished by their ability for long-term self-renewal. Our outcome of
interest is the appearance of a double-mutant cell, which carries a high
potential for leukemic transformation. In our model, a single-hit mutation
carries a slight proliferative advantage over a wild-type stem cells. We
compute marginalized transition probabilities that allow us to capture
important quantitative aspects of our model, including the probability of
observing a double-hit mutant and relevant moments of a single-hit mutation
population over time. We thoroughly explore the model behavior numerically,
varying birth rates across the initial sizes and populations of wild type stem
cells and single-hit mutants, and compare the probability of observing a
double-hit mutant under these conditions. We find that increasing the number
of single-mutants over wild-type particles initially present has a large
effect on the occurrence of a double-mutant, and that it is relatively safe
for single-mutants to be quite proliferative, provided the lentiviral gene
addition avoids creating single mutants in the original insertion process. Our
approach is broadly applicable to an important set of questions in cancer
modeling and other population processes involving multiple stages,
compartments, or types.

###### Keywords:

Mechanistic Models Generating Functions Conditional Moments Hematopoiesis
Stochastic Population Processes Mutagenesis

##  1 Introduction

Branching processes comprise a class of stochastic processes that model how a
collection of discrete individuals changes over time, typically by
reproducing, dying, or transforming type. These processes were first studied
in relation to the extinction of family names among the aristocracy by Francis
Galton and the Reverend Henry Watson kendall1966branching . Branching
processes have since become prevalent across a wide array of scientific
disciplines, from phylogenetics blum2006random  to cancer biology
durrett2015branching  and nuclear physics pazsit2007neutron . In systems
biology contexts related to our interest here, branching processes have been
extensively applied to model clonal expansion in response to tumor suppressor
gene inactivation michor2006stochastic ; nowak2006genetic  and the
evolutionary dynamics of leukemia in response to combination therapy
komarova2005drug ; wodarz2009towards . Stochastic models can be used to study
events such as extinction of a leukemic clone or mutation leading to malignant
transformation.

While many properties of branching processes are well-studied, computing
finite-time quantities that characterize their transient behavior remains
challenging. In this article we propose numerical methods for computing a wide
range of quantities pertinent to a multi-type branching process system
evolving in continuous time. We are motivated by the biological process of
mutation accumulation, specifically in clonal hematopoiesis leading to
myelodysplasia and leukemia, and in gene therapy through lentiviral gene
addition. To this end, we consider a three-type branching process model for
stagewise mutation in a reproducing population of stem cells, and show how a
variety of marginalized transition probabilities that are useful toward
addressing biological questions of interest. This method also enables us to
compute the expected sizes of sub-populations at any time $t$ and related
probabilistic quantities with ease in a general $m$-type model.

A classical mathematical tool used for studying branching processes is the
probability generating function (PGF). The PGF greatly simplifies the
computation of transition probabilities, the fundamental quantities that
determine how a stochastic process behaves. One avenue of attack is to work
with an infinite series of coupled ordinary differential equations (ODEs),
frequently called the Master Equations, that govern the evolution of the
processs states. Instead we exploit the independence of particle lineages to
reduce the computation to solving a small number of ODEs governing the
behavior of the PGFs for processes originating with a single particle of each
type. This property of branching processes—that particles act independently
according to a probability law—makes these simplified mathematical
manipulations possible, greatly simplifying their computational cost.

Our method revolves around three core steps. We first factorize the PGF for
the full process into a product of PGFs of the process beginning with a single
particle of each type. Second, we solve for these PGFs using the pseudo-
generating function technique which involves a simple numerical computation of
a small series of coupled ODEs bailey1990elements  governing the evolution of
the PGF. This is followed by fast numerical inversion based on the Fast
Fourier transform to recover the desired quantities from the PGF solutions.

Much of the prior mathematical literature has largely focused on the “forward”
behavior of branching processes, including the long-run or limiting behavior
of such processes durrett2015branching ; bozic2013evolutionary . The
derivations presented in this article instead help us quantify transient
behavior with an efficient and general numerical recipe, extending these
methods beyond the special cases for which solutions have been previously
established. In particular, finite-time quantities often translate directly to
answer relevant clinical questions, and are necessary within inferential
schemes xu2015likelihood ; xu2019statistical  such as evaluating the marginal
likelihood of data available at a sequence of observation times. Indeed,
transition probabilities fully characterize the probabilistic behavior of
these stochastic models. Our approach targets these and related quantities
applicable to multi-type branching processes beyond our model of mutation
here. This contrasts with previous work that for instance posits a neutral or
symmetric branching process with specific relationships between the birth and
death parameters lan2017fate ; it is worth noting that the more general
version of that model was deemed mathematically intractable, analyzed only in
supplementary simulations. Our contributions fill a methodological gap
suggested in this and related prior literature, which largely resorts to
brute-force simulation and thus requires considerable computation time
lan2017fate .

As a motivating application of these methods, we consider a numerical
examination of gene therapy in hematopoietic stem cell lines. Advances in gene
therapy offer a promising approach to the treatment of genetic hematologic
diseases and immunodeficiency disorders kohn2001gene ; cavazzana2000gene ;
hacein2002sustained . Lentiviral vectors, derived from the human
immunodeficiency virus, can efficiently insert, modify, and delete genes into
a cell’s genome. They have been extensively investigated and optimized over
the last twenty years, and their use is commonplace in academic laboratories
and in gene therapy clinical trials, such as in the treatment of primary
immunodeficiencies and hematologic disorders reinhardt2021long . Most
recently, gene editing through CRISPR technology is being investigated as a
means to treat inherited diseases and offers hope for future translational
work. However, for any gene therapy, when a corrected gene is being inserted
into the genome, it carries a chance of being inserted in a site that is close
to an oncogene. As a result, induction of leukemia remains a major concern for
gene therapy clinical trials kohn2001gene ; marshall2002gene ;
marshall2003second , although more recently developed methods adopt risk
reduction strategies to mitigate these concerns reinhardt2021long . Because
insertion of a corrected gene can disrupt gene structure and alter gene
transcription and gene regulation, there is always a theoretical risk of cell
transformation, although reported events are extremely rare. In the historical
setting of retroviral gene therapy studies, these exceedingly rare events may
occur within weeks of treatment or even up to years later kohn2003occurrence .
Further risk of mutational events in patients undergoing gene therapy clinical
trials comes from the chemotherapy regimens which carry a risk of
hematopoietic stem cell damage and leukemic transformation.

Because leukemia is an unacceptable outcome in any clinical trial, it is
important to quantify the risk of its development in this setting
bonetta2002leukemia . Our goal is to quantify the likelihood of leukemic
transformation and identify conditions under which a leukemic clone is likely
to expand. Using our numerical methods we explore the population dynamics of
transduced hematopoietic stem cells after stem cell transplantation under a
multitype branching process model. This stochastic modeling approach can
quantify the probability of observing rare events under varying clinical
scenarios, such as mutational oncogenesis and engraftment failure that may
occur during stem cell transplant.

We will investigate multiple questions clinically relevant to gene therapy
here, including 1) the time until two mutational events occur leading to
leukemogenesis, and 2) how varying the proliferative advantage of the
transduced population alters the probability of observing leukemic outcomes.
These efforts are motivated by Knudson’s finding that two rate-limiting
mutational steps are required for inactivation of tumor suppressor genes in
sporadic cancers knudsen1971mutation , and apply generally to work studying
initiation and progression of tumors and identifying driver mutations. We
conclude by outlining how our methods may be applied to other clinical
scenarios in cancer research such as the inactivation of the two alleles of a
tumor suppressor gene and the acquisition of multiple driver mutations.

##  2 Mathematical Notation for Branching Processes

We begin with a multi-type branching process whose counts are denoted by the
vector $\mathbf{X}(t)$. A continuous-time branching process $\mathbf{X}(t)$ is
a Markov jump process in which a collection of individuals, in our case cells,
can reproduce and die independently according to a probability distribution.
Each component $X_{i}(t)$ denotes the number of particles of type $i$ at time
$t\geq 0$. Each type has a distinct mean lifespan and reproductive
probabilities, and can give rise to cells of its own type as well as other
types at its time of death.

$X_{1}$$X_{2}$$\lambda_{1}$$\lambda_{2}$$\mu_{1}$$\emptyset$$\mu_{2}$$\emptyset$$\nu$
Figure 1: Reaction diagram for the two-type branching process. $\lambda_{i}$
corresponds to reproduction through binary fission, $\mu_{i}$ corresponds to
removal via cell death, and $\nu$ corresponds to transformation via mutation.

To simplify the exposition and notational load, we focus our attention on the
two-type case depicted in Figure 1 as we present the definitions and
derivations. We then clarify how the methodology applies in the general case,
including the three-type method we implement in our numerical analysis of
mutation in hematopoietic stem cells. Let $a_{i}(k,l)$ denote the rate at
which a particle of type $i$ produces $k$ particles of type one and $l$
particles of type two upon completion of its lifespan. For example,
$a_{1}(2,0)$ is the rate at which a single particle of type one reproduces,
producing two offspring through binary fission. We define the negative total
rate as

$\alpha_{i}\equiv a_{i}(1,0)=-\sum_{(k,l)\neq(1,0)}a_{i}(k,l)$

which satisfies a conservation of total rates $\sum_{(k,l)}a_{i}(k,l)=0$. Our
process is time-homogeneous in that the rates $a_{i}(k,l)$ do not change over
time. This property together with the assumption of particle independence
implies rate-linearity: the reaction rate of any event at time $t$ is given by
the per-particle rate multiplied by the corresponding population at time $t$.

These observations lead us to the instantaneous jump probabilities that
characterize the process. For a short time increment $h$, jumping from a
population of $j$ type one particles to $k$ type one particles and $l$ type
two particles is given by

$\Pr\left(\mathbf{X}(h)=(k,l)\mid\mathbf{X}(0)=(j,0)\right)=j\cdot
a_{1}(k,l)\cdot h+o(h).$ (1)

This form is consistent with our branching process being an instance of a
time-homogeneous continuous time Markov chain lange2010applied . Note that as
$h\downarrow 0$ the instantaneous rate of change is given by $j\cdot
a_{1}(k,l)$ as expected from the previous discusson of rate-linearity. By the
Markov assumption, the lifespan of each particle of type $i$ is exponentially
distributed with rate $-\alpha_{i}$.

For an arbitrary interval of time $t$ that need not be close to zero, we
denote the finite-time transition probabilities

$P_{(i,j),(k,l)}(t)\equiv\Pr\left(\mathbf{X}(t)=(k,l)\mid\mathbf{X}(0)=(i,j)\right).$
(2)

In contrast to (1), the form of $P_{(i,j),(k,l)}(t)$ is typically not readily
available, requiring a large integration step accounting for all possible
process paths and their respective probabilities passing between the endpoints
$\mathbf{X}(0)$ and $\mathbf{X}(t)$ . Towards computing these quantities, we
note that they appear within the probability generating function (PGF), given
by

$\phi_{ij}(t,s_{1},s_{2})\equiv\sum_{k=0}^{\infty}\sum_{j=0}^{\infty}P_{(i,j),(k,l)}(t)s_{1}^{k}s_{2}^{l}=\mathrm{E}\left[s_{1}^{X_{1}(t)}s_{2}^{X_{2}(t)}\mid
X_{1}(0)=i,X_{2}(0)=j\right]$ (3)

where the dummy variables $s_{i}\in\mathbb{C}$ satisfy $\|s_{i}\|\leq 1$ . We
will frequently shorten the notation for the PGF to $\phi_{ij}(t)$,
suppressing dependence on the dummy variables.

The PGF is a primary mathematical object playing a key role in our
computational technique. Notably, given numerical methods to compute the PGF,
it is tractable to invert the PGF using the Fast Fourier Transform to obtain
transition probabilities (2), discussed further below. This is not an entirely
trivial numerical task, however, requiring $O\left(\max(k,l)^{2}\right)$
computations xu2015efficient .

Fortunately, it is often the case that we are concerned primarily with the
distribution of a single sub-population. For instance, in assessing the safety
of a gene therapy we may typically pose questions pertaining to the number of
particles in the mutant or type two population. In such cases, the relevant
series expression can be obtained by marginalizing out all other types from
the joint distribution of the process. Fortunately, doing so becomes
straightforward if one has the joint PGF in hand and further provides the
advantage of reducing the computational complexity considerably. Marginalizing
out the first type, for instance, is accomplished by letting $s_{1}=1$ in the
PGF (3),

$\phi_{ij}(t,s_{1},s_{2})\big{|}_{s_{1}=1}=\sum_{l=0}^{\infty}P_{(i,j),(\cdot,l)}(t)s_{2}^{l}.$

Obtaining these marginalized probabilities is much more computationally
efficient, as it requires inverting only a univariate probability generating
function regardless of the number of total types in the branching process
system. Indeed, though the notation above corresponds to the two-type case,
the analogous quantity pertaining to a type $j$ population in the general
multitype case can be computed by setting $s_{i}=1$ for all populations $i\neq
j$.

##  3 Computing the Probability Generating Function

Here we detail how to compute the probability generating function following a
classical technique detailed in bailey1990elements . We shall generate a
series of Ordinary Differential Equations (ODEs) that govern the evolution of
the PGF over time. We first make the observation that, by particle
independence,

$\phi_{ij}(t)=\phi_{(1,0)}(t)^{i}\phi_{(0,1)}(t)^{j}.$ (4)

For simplicity’s sake we shall shorten the notation
$\phi_{(1,0)}\equiv\phi_{1}$ and $\phi_{(0,1)}\equiv\phi_{2}$. Since each clan
of particles is independent of every other clan, we can view the full PGF as
being composed of the individual contributions from each original ancestor
particle present at time $t=0$. We emphasize here that although the branching
process is defined over a countably infinite state space, by making use of the
master equations applied to the PGF instead of the process itself, we need
only to consider a system of two differential equations. Moreover, this
technique is general and entails $m$ differential equations for the $m$-type
case, which we discuss after detailing the two-type case in order to lighten
notation.

Equation (4) suggests it suffices to work with the PGFs of the process
starting with a particle of each given type. Again illustrating with two-type
notation, consider the small-time expansion of $\phi_{1}$ using (1),

$\displaystyle\phi_{1}$
$\displaystyle(t,s_{1},s_{2})=\mathrm{E}\left[s_{1}^{X_{1}(t)}s_{2}^{X_{2}(t)}\mid
X_{1}(0)=1,X_{2}(0)=0\right]=\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}P_{(1,0),(k,l)}(t)s_{1}^{k}s_{2}^{l}$
$\displaystyle=\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}\left[\mathbf{1}_{k=1,l=0}+a_{1}(k,l)t+o(t)\right]s_{1}^{k}s_{2}^{l}=s_{1}+t\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}a_{1}(k,l)s_{1}^{k}s_{2}^{l}+o(t)$
(5)

as $t\downarrow 0$. Here $\mathbf{1}$ denotes the indicator function. Above,
we see that the final expression in Equation (3) involves the pseudo-
generating function bailey1990elements

$u_{i}(s_{1},s_{2})=\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}a_{i}(k,l)s_{1}^{k}s_{2}^{l}.$
(6)

Using classical Chapman-Kolmogorov arguments with details in the Appendix, we
arrive at the system of backward equations

$\displaystyle\frac{\ \text{d}missing}{\ \text{d}t}\phi_{1}(t)$
$\displaystyle=u_{1}(\phi_{1}(t),\phi_{2}(t)),\qquad\phi_{1}(0)=s_{1}$
$\displaystyle\frac{\ \text{d}missing}{\ \text{d}t}\phi_{2}(t)$
$\displaystyle=u_{2}(\phi_{1}(t),\phi_{2}(t)),\qquad\phi_{2}(0)=s_{2}.$

The initial conditions follow from the definition of the PGF (3).

We will now focus on the specific rates $a_{i}(k,l)$ that are used in the two-
type branching process, listed below:

$\displaystyle a_{1}(2,0)=\lambda_{1},\qquad a_{1}(0,0)=\mu_{1},\qquad
a_{1}(0,1)=\nu,\qquad a_{1}(1,0)=-(\lambda_{1}+\mu_{1}+\nu)$ $\displaystyle
a_{2}(0,2)=\lambda_{2},\qquad a_{2}(0,0)=\mu_{2},\qquad
a_{2}(0,1)=-(\lambda_{2}+\mu_{2}).$

Writing these in terms of pseudo-generating functions, we have

$\displaystyle u_{1}(s_{1},s_{2})=\lambda_{1}s_{1}^{2}+\nu
s_{2}-(\lambda_{1}+\mu_{1}+\nu)s_{1}+\mu_{1}$ $\displaystyle
u_{2}(s_{1},s_{2})=\lambda_{2}s_{2}^{2}-(\lambda_{2}+\mu_{2})s_{2}+\mu_{2}.$

We use these expressions in the differential equations for the single-particle
generating functions to arrive at the desired backwards equations and initial
conditions, as follows:

$\displaystyle\frac{\ \text{d}missing}{\ \text{d}t}\phi_{1}(t)$
$\displaystyle=\lambda_{1}\phi_{1}^{2}+\nu\phi_{2}-(\lambda_{1}+\mu_{1}+\nu)\phi_{1}+\mu_{1},\qquad$
$\displaystyle\phi_{1}(0)=s_{1}$ (7) $\displaystyle\frac{\ \text{d}missing}{\
\text{d}t}\phi_{2}(t)$
$\displaystyle=\lambda_{2}\phi_{2}^{2}-(\lambda_{2}+\mu_{2})\phi_{2}+\mu_{2},\qquad\qquad$
$\displaystyle\phi_{2}(0)=s_{2}.$ (8)

###  3.1 Solving for $\phi_{2}$

Upon inspection, the differential equation (8) is a nonlinear first order
Riccati equation, which is fortunately solvable. This is because the ansatz
$\phi_{2}=K$ provides us the constant particular solutions $\phi_{2}=1$ and
$\phi_{2}=\mu_{2}/\lambda_{2}$. As is standard for solving such a Riccati
equation, we make the substitution $z=\frac{1}{\phi_{2}-1}$, implying that
$\phi_{2}=1+1/z$. Making this substitution in (8) gives

$\displaystyle\phi_{2}^{\prime}=-\frac{z^{\prime}}{z^{2}}$
$\displaystyle=\mu_{2}-(\lambda_{2}+\mu_{2})\left(1+\frac{1}{z}\right)+\lambda_{2}\left(1+\frac{1}{z}\right)^{2}$
$\displaystyle=\mu_{2}-(\lambda_{2}+\mu_{2})-\frac{\lambda_{2}+\mu_{2}}{z}+\lambda_{2}\left(\frac{1}{z^{2}}+\frac{2}{z}+1\right)$
$\displaystyle=\frac{\lambda_{2}-\mu_{2}}{z}+\frac{\lambda_{2}}{z^{2}}.$

Multiplying through by $-z^{2}$ and rearranging gives a linear first order
differential that is straightforward to solve using an integrating factor,

$\displaystyle z^{\prime}$ $\displaystyle=(\mu_{2}-\lambda_{2})z-\lambda_{2}$
$\displaystyle z$
$\displaystyle=\frac{\lambda_{2}}{\mu_{2}-\lambda_{2}}+Ce^{(\mu_{2}-\lambda_{2})t}$

for some constant $C$ that depends on the initial conditions. We substitute
$z$ back into $\phi_{2}=1+1/z$ and include the initial conditions
$\phi_{2}(0)=s_{2}$ to get the solution

$\phi_{2}(t,s_{2})=1+\left[\frac{\lambda_{2}}{\mu_{2}-\lambda_{2}}+\left(\frac{1}{s_{2}-1}+\frac{\lambda_{2}}{\lambda_{2}-\mu_{2}}\right)e^{(\mu_{2}-\lambda_{2})t}\right]^{-1}.$
(9)

Note that in the critical case $\lambda_{2}=\mu_{2}$, Equation (9) becomes

$\phi_{2}(t,s_{2})=1+\left[\frac{1}{s_{2}-1}-\lambda_{2}t\right]^{-1}.$

###  3.2 Solving for $\phi_{1}$ and the general case

The differential equation (7) governing $\phi_{1}$ is more difficult due to
the presence of a mutation term involving $\phi_{2}$. We now have a single
inhomogeneous nonlinear Riccati equation that cannot be solved by a simple
substitution trick. However, we note that it is straightforward to solve
$\phi_{1}$ numerically using any number of differential equation solvers; we
use `DifferentialEquations.jl` package in the `Julia` programming language
rackauckas2017differentialequations ; bezanson2017julia . The computational
bottleneck to evaluate the joint PGF at any time $t>0$ amounts to numerically
solving a single ODE.

$X_{1}$$X_{2}$$X_{3}$$\lambda_{1}$$\lambda_{2}$$\lambda_{3}$$\mu_{1}$$\emptyset$$\mu_{2}$$\emptyset$$\mu_{3}$$\emptyset$$\nu_{1}$$\nu_{2}$
Figure 2: Reaction diagram for the three-type branching process. $\lambda_{i}$
corresponds to reproduction through binary fission, $\mu_{i}$ denotes removal
via death, and $\nu_{i}$ represents transformation via mutation.

We have focused on the two-type model to provide a clear exposition to bridge
classical methods to the numerical techniques we will utilize, but these
immediately extend beyond the two-type case. Given an instance of an $m$-type
model in the analogous framework, it is straightforward to write the system of
backward equations analogous to (7) and (8). Typically, the equation
corresponding to the final type may be solvable, while we recommend
numerically solving the remaining system of differential equations. For
example, the survival probability of a birth-death model has an analytical
solution.

Concretely, we illustrate how this approach readily applies to a three-stage
model that will be used to study mutation in a hematopoietic model in Section
4. Figure 2 illustrates the reaction rates characterizing the model. Mirroring
the notation earlier, we abbreviate $\phi_{1}=\phi_{(1,0,0)}$,
$\phi_{2}=\phi_{(0,1,0)}$, and $\phi_{3}=\phi_{(0,0,1)}$. The differential
equations that govern the single-particle generating functions are thus:

$\displaystyle\frac{d}{dt}\phi_{1}$
$\displaystyle=\lambda_{1}\phi_{1}^{2}-\left(\lambda_{1}-\mu_{1}-\nu_{1}\right)\phi_{1}+\nu_{1}\phi_{2}+\mu_{1}$
$\displaystyle\frac{d}{dt}\phi_{2}$
$\displaystyle=\lambda_{2}\phi_{2}^{2}-\left(\lambda_{2}-\mu_{2}-\nu_{2}\right)\phi_{2}+\nu_{2}\phi_{3}+\mu_{2}$
$\displaystyle\frac{d}{dt}\phi_{3}$
$\displaystyle=\lambda_{3}\phi_{3}^{2}-\left(\lambda_{3}-\mu_{3}\right)\phi_{3}+\mu_{3}.$

Note that $\frac{d}{dt}\phi_{3}$ has the same structure as
$\frac{d}{dt}\phi_{2}$ from the two-type exposition (8); we can use the same
argument to show that

$\phi_{3}(t,s_{3})=1+\left[\frac{\lambda_{3}}{\mu_{3}-\lambda_{3}}+\left(\frac{1}{s_{3}-1}+\frac{\lambda_{3}}{\lambda_{3}-\mu_{3}}\right)e^{(\mu_{3}-\lambda_{3})t}\right]^{-1},$

and likewise solve $\phi_{1}$ and $\phi_{2}$ numerically.

It is worth noting that in special cases such as two-type models described
above, it is possible to derive exact solutions antal2011exact . From a
numerical perspective, terms such as hypergeometric functions arising in these
such solutions can present numerical challenges in some parameter regimes
(xu2015likelihood, , Appendix A). When this is the case, which often occurs
when these solutions are required within iterative schemes, it can be
preferable to solve $\phi_{1}$ using numerical differential equation solvers.
Moreover, the numerical recipe we advocate above for the solving the system
holds for any $m$-type model in our class.

###  3.3 Inverting the Numerical PGF

Once we have solved for $\phi_{1}$ and $\phi_{2}$ to a specified time point,
we can compute the full PGF $\phi_{ij}$ via (4). We then use a technique
introduced in lange1982calculation  to invert the PGF to obtain our desired
marginalized univariate probability function. We begin with a change of
variables $s=e^{2\pi iw}$, placing our PGF argument $s$ on the unit circle in
the complex plane. This transforms the generating function into a periodic
function,

$\phi_{ij}\left(t,e^{2\pi iw}\right)=\sum_{k=0}^{\infty}c_{k}(t)e^{2\pi iwk}$

where $c_{k}(t)$ are the coefficients that contain our desired transition
probabilities and $i=\sqrt{-1}$. $c_{k}$ is also the $k$th coefficient of a
Fourier series, and therefore can be obtained by a straightforward inversion

$c_{k}(t)=\int_{0}^{1}\phi_{ij}\left(t,e^{2\pi iw}\right)e^{-2\pi iwk}dw.$

In practice, the integral can be evaluated numerically via discretization. A
Riemann sum approximation is appropriate for large $N$, yielding

$c_{k}(t)\equiv
P_{(ij),k}(t)\approx\frac{1}{N-1}\sum_{l=0}^{N-1}\phi_{ij}\left(t,e^{2\pi
il/N}\right)e^{-2\pi ilk/N}.$

Note that one can increase the accuracy of the approximation by increasing the
size of the gridlength, i.e. make $N$ larger. Most importantly, rather than
naively summing all terms in the Riemann sum, one can apply the Fast Fourier
Transform (FFT) to recover all of the coefficients for $k=0$ to $N-1$
simultaneously in a single computation. This yields the marginalized
transition probabilities in a numerically stable and efficient manner
regardless of the number of types. In contrast, without using the marginalized
forms, this technique would involve summations over each population type,
dramatically increasing the computational cost.

###  3.4 Extracting Relevant Probabilities and Expectations

With a numerical solution for $\phi_{ij}$ now readily available, obtaining
expectations becomes straightforward via numerical differentiation. For
instance, from starting population $(m,n)$, one can compute the average type
$i$ population after any time interval of length $t$ as

$\mathrm{E}\left[X_{i}(t)\mid\mathbf{X}(0)=(m,n)\right]=\frac{\partial}{\partial
s_{i}}\phi_{(m,n)}(t,s_{1},s_{2})\big{|}_{s_{1}=s_{2}=1}.$

Higher order moments can be obtained in a similar fashion by computing further
partial derivatives; see lange2010applied .

While such quantities that are expressible as expectations can also be
computed by way of moment generating functions, other more nuanced quantities
that translate to specific, interpretable questions often require the full
PGF. Many clinical questions can be stated in terms of thresholds, such as the
probability that a double-mutant population exceeds a critical level. These
probabilities can be obtained by summing on the probability generating
function of the population of interest after the threshold, while
marginalizing over all other subtypes:

$\Pr\left(X_{2}(t)\geq
M\mid\mathbf{X}(0)=(i,j)\right)=\sum_{k=0}^{\infty}\sum_{l=M}^{\infty}P_{(i,j),(k,l)}(t).$

Similarly, for populations with more than two types, marginalizing over
multiple types or species $i$ is accomplished by substituting $s_{i}=1$ into
the joint PGF.

Indeed, when the question of interest regards a critical level pertaining to
one type, we wish to marginalize out all but one of the populations.
Specifically marginalizing over $X_{1}$ can be done via

$\phi_{ij}(t,s_{1},s_{2})\big{|}_{s_{1}=1}=\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}P_{(i,j),(k,l)}(t)s_{2}^{l}=\sum_{l=0}^{\infty}P_{(i,j),(\cdot,l)}(t)s_{2}^{l}.$

This marginalized PGF now describes the probabilistic behavior of the $X_{2}$
population, having correctly accounted for the probabilistic behavior of all
configurations of the rest of the system. As such it becomes a univariate
function in $s_{2}$. The previous numerical FFT approach to invert the
generating function becomes computationally efficient to obtain the marginal
transition probabilities $P_{(i,j),(\cdot,l)}(t)$ for $X_{2}(t)$. Then we
compute the desired threshold probabilities as

$\Pr\left(X_{2}(t)\geq
M\mid\mathbf{X}(0)=(i,j)\right)=1-\sum_{l=0}^{M-1}P_{(i,j),(\cdot,l)}(t).$

##  4 Numerical Results

We now investigate the three-type branching process as a model for diploid
mutations in a growing population of hematopoietic stem cells. This model can
be applied to a growing population of transplanted stem cells in gene therapy,
particularly to the acquisition of mutations that lead to oncogenesis. The
rate parameters are listed in Table 1, corresponding to the reactions depicted
previously in Figure 2.

Rate Symbol |  | Diagram |  | Parameter Value |
---|---|---|---|---|---
$\lambda_{1}$ |  | $X_{1}\rightarrow X_{1}+X_{1}$ |  | $2.4\times 10^{-2}$ |
$\mu_{1}$ |  | $X_{1}\rightarrow\emptyset$ |  | $1.4\times 10^{-2}$ |
$\nu_{1}$ |  | $X_{1}\rightarrow X_{2}$ |  | $\lambda_{1}\times 10^{-8}$ |
$\lambda_{2}$ |  | $X_{2}\rightarrow X_{2}+X_{2}$ |  | Varies |
$\mu_{2}$ |  | $X_{2}\rightarrow\emptyset$ |  | $1.4\times 10^{-2}$ |
$\nu_{2}$ |  | $X_{2}\rightarrow X_{3}$ |  | $\lambda_{2}\times 10^{-8}$ |
$\lambda_{3}$ |  | $X_{3}\rightarrow X_{3}+X_{3}$ |  | Varies |
$\mu_{3}$ |  | $X_{3}\rightarrow\emptyset$ |  | $1.4\times 10^{-2}$ |
Table 1: Reactions and parameters for the three-type branching process. The
rate parameters correspond to units of 1/week and are sourced from
abkowitz2002evidence ; pillay2010vivo ; patel2017fate ; michor2005dynamics .

We begin by comparing how varying the initial counts of the wild-type and
single-mutants affects the probability of observing at least one double-mutant
one year after transplant. Figure 3 depicts these probabilities and draws the
expected conclusion that increasing the number of single-mutants over wild-
type particles initially present has a larger effect on the occurence of a
double-mutant. This result has clinical relevance as it places a bound on how
mutagenic the lentiviral gene addition can be before deleterious effects can
be expected. If the mutagenesis is rare, producing only one to ten single-
mutants in the transplanted population, then there is a very small risk that a
deleterious second mutation event will occur within a year. Additionally, the
number of transplanted wild-type stem cells can be quite large before
mutational events become likely, which bodes well for transplanting large
populations of stem cells.

Figure 3: Log probability of observing at least one double-mutant at $t=52$
weeks while varying the initial counts of $X_{1}$ and $X_{2}$.

Next we compare varying the birth rates across the population of wild-type
stem cells and single-mutants. Figure 4 assumes no initial single-mutants, and
so the birth rate of the wild-type population $X_{1}$ is rate-limiting on the
probability of observing a double-mutant. We observe that changing the wild-
type birth rate by a factor of ten substantially increases the odds of
observing a double-mutant, while changing the single-mutant birth rate has a
much smaller effect in comparison. This finding suggests it is relatively safe
for the single-mutations to be quite proliferative if the lentiviral gene
addition avoids creating any single-mutants in the original insertion process.

Figure 4: Log probability of observing at least one double-mutant at $t=52$
weeks while varying the birth rates $\lambda_{i}$ of $X_{1}$ and $X_{2}$. The
probability is depicted in log scale. We have $\lambda_{i}=0.024\times
10^{y_{i}}$, e.g. when $y_{i}=0$, $\lambda_{i}=0.024$. The initial counts are
$X_{1}(0)=1000$ and $X_{2}(0)=0$.

Since we are solving the ODEs for the univariate PGFs forward in time, we can
efficiently examine how the probability of observing a double-mutant changes
over time based off of varying initial population counts, without any
additional computation. This information is immediately available from the
intermediate steps of the ODE solver, visualized in Figure 5. Qualitatively,
we see three phases of growth corresponding to a fast initial phase, followed
by an exponential growth in the probability, and finally tapering out towards
one as a double-mutant becomes guaranteed to be observed. We also observe that
the duration of the second phase depends on the initial population count.

Figure 5: Probability of observing at least one double-mutant over time while
varying the initial counts of $X_{1}$ and $X_{2}$. The probability is depicted
in log scale.

##  5 Conclusions

This paper proposes ways to numerically evaluate marginal probabilistic
quantities based on the PGF $\phi_{ijk}$ of a three-type branching process.
Our technique is exact up to numerical error in inverting the FFT integral,
and we show how many clinically relevant questions can be reduced to tractable
univariate problems. For instance, in our experiments we are concerned
primarily with the double-mutant population $X_{3}$, and our method yields
information about the probabilities of leukemogenesis after gene therapy. This
presents important information for clinicians about how leukemia may develop
under a variety of conditions including mutagenesis from the initial gene
insertion, as well as through differing proliferative rates among the mutant
and wild-type populations.

Our methodology is general in that we can alter the structure of the
underlying mathematical process without significantly changing the numerical
processes used to obtain the transition probabilities. For example, both the
shift process $X_{1}\rightarrow X_{2}$ and the mutational birth process
$X_{1}\rightarrow X_{1}+X_{2}$ can be incorporated into our system of
equations. These changes will change the structure of the underlying
differential equations to be solved, but they do not change the overall
numerical manipulations required to solve for and invert the probability
generating function. This direct numerical technique provides an attractive
alternative to computationally intensive simulation-based approaches. This
computational efficiency becomes critical in studies requiring exploration of
many parameters, iterative optimization, or within statistical inferential
schemes.

Our approach applies to any multitype branching process model under the
standard rate-linearity assumption, enabling us to calculate a wide variety of
pertinent transition probabilities, expectations and higher moments toward
exploring the transient behavior of the particles. We see that the framework
can be used to answer questions about the population of any given type
$X_{i}(t)$ in a general multitype process. Similarly, one can pose related
questions related to sums or differences of populations, such as the
probability that one population becomes dominant or outnumbers another
(lange2010applied, , see 13.3.1). In some cases, we can establish connections
to related models and results in the literature. For instance, an anonymous
reviewer remarks that we may recognize the equation corresponding to the final
type as the birth-death equations, for which additional tools are available.
Indeed, if one investigates the probability of a double-mutant in the three
stage model depicted in Figure 2, the problem can be equivalently case in
terms of survival probabilities of birth-death processes can be used
denes1996exact ; crawford2012transition . Exploring further connections and
possible extensions of these techniques to systems with non-linear interaction
dynamics will be a fruitful avenue for future work.

## Appendix

To derive the backward equations, we begin by rearranging (3), divide by $t$,
then send $t\downarrow 0$ to get an expression for the derivatives of the
single-particle PGFs in terms of pseudo-generating functions (6),

$\frac{\ \text{d}\phi_{1}(t,s_{1},s_{2})}{\
\text{d}t}\Big{|}_{t=0}=u_{1}(s_{1},s_{2}),\qquad\frac{\
\text{d}\phi_{2}(t,s_{1},s_{2})}{\
\text{d}t}\Big{|}_{t=0}=u_{2}(s_{1},s_{2}).$ (10)

Now we may apply the Chapman-Kolmogorov equations for $\phi_{1}$ along with
(4) to find that

$\displaystyle\phi_{1}(t+h,s_{1},s_{2})$
$\displaystyle=\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}P_{(1,0),(k,l)}(t+h)s_{1}^{k}s_{2}^{l}$
$\displaystyle=\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}\left[\sum_{i=0}^{\infty}\sum_{j=0}^{\infty}P_{(1,0),(i,j)}(t)P_{(i,j),(k,l)}(h)\right]s_{1}^{k}s_{2}^{l}$
$\displaystyle=\sum_{i=0}^{\infty}\sum_{j=0}^{\infty}P_{(1,0),(i,j)}(t)\left[\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}P_{(i,j),(k,l)}(h)s_{1}^{k}s_{2}^{l}\right]$
$\displaystyle=\sum_{i=0}^{\infty}\sum_{j=0}^{\infty}P_{(1,0),(i,j)}(t)\phi_{ij}(h,s_{1},s_{2})$
$\displaystyle=\sum_{i=0}^{\infty}\sum_{j=0}^{\infty}P_{(1,0),(i,j)}(t)\phi_{1}(h,s_{1},s_{2})^{i}\phi_{2}(h,s_{1},s_{2})^{j}$
$\displaystyle=\phi_{1}\big{(}t,\phi_{1}(h,s_{1},s_{2}),\phi_{2}(h,s_{1},s_{2})\big{)}.$

The relationship is time-symmetric, therefore we also have

$\phi_{1}(t+h,s_{1},s_{2})=\phi_{1}\big{(}h,\phi_{1}(t,s_{1},s_{2}),\phi_{2}(t,s_{1},s_{2})\big{)}.$
(11)

Naturally the same relationships exist for $\phi_{2}$.

We now derive the backward Kolmogorov equations for $\phi_{1}$ and $\phi_{2}$,
obtained from a Taylor expansion of $\phi_{1}$ around $t$ together with (10)
and (11),

$\displaystyle\phi_{1}(t+h,s_{1},s_{2})$
$\displaystyle=\phi_{1}(t,s_{1},s_{2})+\frac{\
\text{d}\phi_{1}(t+h,s_{1},s_{2})}{\ \text{d}t}\Big{|}_{h=0}h+o(h)$
$\displaystyle=\phi_{1}(t,s_{1},s_{2})+\frac{\
\text{d}\phi_{1}\big{(}h,\phi_{1}(t,s_{1},s_{2}),\phi_{2}(t,s_{1},s_{2})\big{)}}{\
\text{d}t}\Big{|}_{h=0}h+o(h)$
$\displaystyle=\phi_{1}(t,s_{1},s_{2})+u_{1}\big{(}\phi_{1}(t,s_{1},s_{2}),\phi_{2}(t,s_{1},s_{2})\big{)}h+o(h).$

The same actions produce an analogous expression for $\phi_{2}$. Performing
the previous rearrangement, dividing by $h$, and sending $h$ to zero gives us
our desired system of backwards Kolmogorov equations.

## References

  * (1) Abkowitz, J.L., Catlin, S.N., McCallie, M.T., Guttorp, P.: Evidence that the number of hematopoietic stem cells per animal is conserved in mammals.  Blood, The Journal of the American Society of Hematology 100(7), 2665–2667 (2002)
  * (2) Antal, T., Krapivsky, P.: Exact solution of a two-type branching process: models of tumor progression.  Journal of Statistical Mechanics: Theory and Experiment 2011(08), P08018 (2011)
  * (3) Bailey, N.T.: The elements of stochastic processes with applications to the natural sciences, vol. 25.  John Wiley & Sons (1990)
  * (4) Bezanson, J., Edelman, A., Karpinski, S., Shah, V.B.: Julia: A fresh approach to numerical computing.  SIAM Review 59(1), 65–98 (2017)
  * (5) Blum, M.G., François, O.: Which random processes describe the tree of life? A large-scale study of phylogenetic tree imbalance.  Systematic Biology 55(4), 685–691 (2006)
  * (6) Bonetta, L.: Leukemia case triggers tighter gene-therapy controls.  Nature Medicine 8(11), 1189–1189 (2002)
  * (7) Bozic, I., Reiter, J.G., Allen, B., Antal, T., Chatterjee, K., Shah, P., Moon, Y.S., Yaqubie, A., Kelly, N., Le, D.T., et al.: Evolutionary dynamics of cancer in response to targeted combination therapy.  elife 2, e00747 (2013)
  * (8) Cavazzana-Calvo, M., Hacein-Bey, S., de Saint Basile, G., Gross, F., Yvon, E., Nusbaum, P., Selz, F., Hue, C., Certain, S., Casanova, J.L., et al.: Gene therapy of human severe combined immunodeficiency (scid)-x1 disease.  Science 288(5466), 669–672 (2000)
  * (9) Crawford, F.W., Suchard, M.A.: Transition probabilities for general birth–death processes with applications in ecology, genetics, and evolution.  Journal of mathematical biology 65(3), 553–580 (2012)
  * (10) Denes, J., Krewski, D.: An exact representation for the generating function for the moolgavkar-venzon-knudson two-stage model of carcinogenesis with stochastic stem cell growth.  Mathematical biosciences 131(2), 185–204 (1996)
  * (11) Durrett, R.: Branching process models of cancer.  In: Branching process models of cancer, pp. 1–63. Springer (2015)
  * (12) Hacein-Bey-Abina, S., Le Deist, F., Carlier, F., Bouneaud, C., Hue, C., De Villartay, J.P., Thrasher, A.J., Wulffraat, N., Sorensen, R., Dupuis-Girod, S., et al.: Sustained correction of x-linked severe combined immunodeficiency by ex vivo gene therapy.  New England Journal of Medicine 346(16), 1185–1193 (2002)
  * (13) Kendall, D.G.: Branching processes since 1873.  Journal of the London Mathematical Society 1(1), 385–406 (1966)
  * (14) Knudsen, A.G.J.: Mutation and cancer: Statistical study of retinoblastoma.  Proceedings of the National Academy of Sciences of the United States of America 68(4), 820–823 (1971)
  * (15) Kohn, D.B.: Gene therapy for genetic haematological disorders and immunodeficiencies.  Journal of internal medicine 249(4), 379–390 (2001)
  * (16) Kohn, D.B., Sadelain, M., Glorioso, J.C.: Occurrence of leukaemia following gene therapy of x-linked scid.  Nature Reviews Cancer 3(7), 477–488 (2003)
  * (17) Komarova, N.L., Wodarz, D.: Drug resistance in cancer: principles of emergence and prevention.  Proceedings of the National Academy of Sciences 102(27), 9714–9719 (2005)
  * (18) Lan, X., Jörg, D.J., Cavalli, F.M., Richards, L.M., Nguyen, L.V., Vanner, R.J., Guilhamon, P., Lee, L., Kushida, M.M., Pellacani, D., et al.: Fate mapping of human glioblastoma reveals an invariant stem cell hierarchy.  Nature 549(7671), 227–232 (2017)
  * (19) Lange, K.: Calculation of the equilibrium distribution for a deleterious gene by the finite fourier transform.  Biometrics pp. 79–86 (1982)
  * (20) Lange, K.: Applied probability.  Springer Science & Business Media (2010)
  * (21) Marshall, E.: Gene therapy a suspect in leukemia-like disease.  Science 298(5591), 34–35 (2002)
  * (22) Marshall, E.: Second child in french trial is found to have leukemia.  Science 299(5605), 320–320 (2003)
  * (23) Michor, F., Hughes, T.P., Iwasa, Y., Branford, S., Shah, N.P., Sawyers, C.L., Nowak, M.A.: Dynamics of chronic myeloid leukaemia.  Nature 435(7046), 1267–1270 (2005)
  * (24) Michor, F., Nowak, M.A., Iwasa, Y.: Stochastic dynamics of metastasis formation.  Journal of Theoretical Biology 240(4), 521–530 (2006)
  * (25) Nowak, M.A., Michor, F., Iwasa, Y.: Genetic instability and clonal expansion.  Journal of theoretical biology 241(1), 26–32 (2006)
  * (26) Patel, A.A., Zhang, Y., Fullerton, J.N., Boelen, L., Rongvaux, A., Maini, A.A., Bigley, V., Flavell, R.A., Gilroy, D.W., Asquith, B., et al.: The fate and lifespan of human monocyte subsets in steady state and systemic inflammation.  Journal of Experimental Medicine 214(7), 1913–1923 (2017)
  * (27) Pázsit, I., Pál, L.: Neutron fluctuations: A treatise on the physics of branching processes.  Elsevier (2007)
  * (28) Pillay, J., Den Braber, I., Vrisekoop, N., Kwast, L.M., De Boer, R.J., Borghans, J.A., Tesselaar, K., Koenderman, L.: In vivo labeling with 2h2o reveals a human neutrophil lifespan of 5.4 days.  Blood, The Journal of the American Society of Hematology 116(4), 625–627 (2010)
  * (29) Rackauckas, C., Nie, Q.: Differentialequations. jl–a performant and feature-rich ecosystem for solving differential equations in julia.  Journal of Open Research Software 5(1) (2017)
  * (30) Reinhardt, B., Habib, O., Shaw, K.L., Garabedian, E., Carbonaro-Sarracino, D.A., Terrazas, D., Fernandez, B.C., De Oliveira, S., Moore, T.B., Ikeda, A.K., et al.: Long-term outcomes after gene therapy for adenosine deaminase severe combined immune deficiency (ada scid).  Blood (2021)
  * (31) Wodarz, D., Komarova, N.: Towards predictive computational models of oncolytic virus therapy: basis for experimental validation and model selection.  PloS one 4(1), e4271 (2009)
  * (32) Xu, J., Guttorp, P., Kato-Maeda, M., Minin, V.N.: Likelihood-based inference for discretely observed birth–death-shift processes, with applications to evolution of mobile genetic elements.  Biometrics 71(4), 1009–1021 (2015)
  * (33) Xu, J., Koelle, S., Guttorp, P., Wu, C., Dunbar, C., Abkowitz, J.L., Minin, V.N., et al.: Statistical inference for partially observed branching processes with application to cell lineage tracking of in vivo hematopoiesis.  The Annals of Applied Statistics 13(4), 2091–2119 (2019)
  * (34) Xu, J., Minin, V.N.: Efficient transition probability computation for continuous-time branching processes via compressed sensing.  In: Uncertainty in artificial intelligence: proceedings of the… conference. Conference on Uncertainty in Artificial Intelligence, vol. 2015, p. 952. NIH Public Access (2015)