Numerical approximations of population size distributions for multi-type branching processes

Abstract

Continuous-time multi-type branching processes are fundamental models for expanding and migrating populations with cancer evolution being a prototypical example. Inferring model parameters, like mutation and growth rates, from time-series count data requires efficient computation of population size distributions. Existing methods are mainly based on large-time or large-number asymptotics, which rely on either restricted initial conditions or simplified interactions between cell types. Here, we introduce two numerical approximations of population size distributions for multi-type branching processes on directed graphs with arbitrary initialization. The first approach combines a saddle-point approximation with numerical integration of the probability generating function. We characterize admissibility and establish conditions for saddle-point existence and uniqueness. For directed acyclic graphs, the second approach provides a large-time small-mutation-rate alternative based on closed-form approximate Laplace transforms and efficient numerical inversion. We benchmark the accuracy and speed of both solutions in simulations, showing substantial improvement over the state-of-the-art large-number approximation and orders of magnitude speedup over Gillespie's stochastic simulation algorithm at matching accuracy. We apply our methods to analyze the relapse dynamics of an acute myeloid leukemia patient, where rapid parameter scans over a six-type patient-specific mutation tree quantify how unobserved remission burden and treatment-altered fitness can explain relapse. Our methods provide computational building blocks for future likelihood-based inference in cancer evolution and other expanding populations.