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

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.

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.

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

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

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

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

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),

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$.

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,

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),

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

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

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:

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

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:

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.

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.

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.

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.

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.

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),

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

The relationship is time-symmetric, therefore we also have

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),

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.