My current research interests fall into several main areas (click on the links to see descriptions of some of my projects in these areas):
- Mathematical models of cancer evolution
- Evolution of drug resistance in cancer and optimal treatment strategies
- Numerical methods for stochastic partial differential equations

Mathematical models of cancer evolution

Cancer can be described as a genetic evolutionary process within the body. During its progression, a clone arising from a single cancer-initiating cell with misregulated growth or death signals multiplies exponentially to become a neoplasm consisting of 10^12 cells or greater. I am interested in stochastic modeling of phenomena such as cancer initiation from healthy tissue and its dependence on tissue structure, the establishment of genetic heterogeneity within a tumor, and the dynamics of drug response. I am also interested in using dynamical models to study the effects of a cancer cell differentiation hierarchy on tumorigenesis and treatment outcomes in various cancer types.

Eradication of Chronic Myeloid Leukemia Stem Cells. Imatinib mesylate (Gleevec) is currently the standard treatment for chronic myeloid leukemia (CML) and elicits a large reduction in leukemic cell burden in most patients. However, strong evidence suggests that imatinib does not cure the disease; approximately 20% of patients relapse within three years, and discontinuation of imatinib therapy often leads to a rebound of the leukemic cell burden. It has been postulated that the disease outcome may be improved by administering imatinib in conjunction with the Granulocyte-Colony Stimulating Factor (G-CSF), a growth factor which wakes up the quiescent stem cells and sensitizes them to imatinib. In this study, we designed a novel mathematical model of stem cell quiescence to investigate the treatment response to imatinib and G-CSF. We found that adding G-CSF to an imatinib treatment protocol leads to observable effects only if the majority of leukemic stem cells are quiescent. Our model also predicted that adding G-CSF leads to a higher risk of resistance, since it increases the number of leukemic stem cell divisions and thus the probability of acquiring a resistance mutation. PDF (w/ M. Drummond, B. Clarkson, T. Holyoake, F. Michor)

Evolutionary dynamics of tumorigenesis with random mutational fitness advances. After a cancer cell initiates clonal expansion, the processes of mutation and selection occur within an exponentially growing population. In this work we model this clonal expansion phase of cancer progression with a multi-type branching process. Each time a cell divides, a mutation may arise with small probability which produces a new cell type, conferring a random increment to the daughter cell's fitness. Under this model, we analyzed the asymptotic rate of expansion for the k-th generation of mutants, for both bounded and unbounded mutational fitness distributions. PDF (w/R. Durrett, K. Leder, J. Mayberry, F. Michor)


Evolution of intra-tumor heterogeneity.
To study the evolution of heterogeneity during tumorigenesis, we investigated a mathematical model of tumorigenesis using multi-type branching processes with random mutational fitness advances. We considered several ecological measures of diversity in the population, such as the Simpson's Index and the contribution of the largest clones to the population. We characterized the Simpson's Index for the first wave of mutants in the population, and in particular, its dependence on the parameters of the mutational fitness distribution. We also investigated the heterogeneity between successive waves of mutants and the arrival times of each wave in the small mutation limit. PDF (w/R. Durrett, K. Leder, J. Mayberry, F. Michor)

Stochastic dynamics of cancer initiation. Elucidation of the dynamics of cancer initiation is of importance for an understanding of tumor evolution and cancer incidence data. Here we developed a mathematical framework to study the processes of cancer initiation. Cells at risk of accumulating oncogenic mutations are organized into small compartments of cells and proliferate according to a stochastic process. During each cell division, an (epi)genetic alteration may arise which leads to a random fitness change, drawn from a probability distribution. Cancer is initiated when a cell gains a fitness sufficiently high to escape from the homeostatic mechanisms of the cell compartment. To investigate cancer initiation during a human lifetime, we considered a `race' between this fitness process and the aging process of the person; the latter is modeled as a second stochastic Markov process in an aging dimension. This model allows us to investigate the dynamics of cancer initiation and its dependence on the mutational fitness distribution. This framework also provides a methodology to assess the effects of different life expectancy distributions on lifetime cancer incidence. PDF (w/K. Leder, F. Michor)

Evolution of drug resistance and optimal treatment strategies

Anti-cancer drugs targeted to specific oncogenic pathways have shown promising therapeutic results in the past few years; however, drug resistance remains an important obstacle for these therapies as well as for conventional cytotoxic drugs (such as chemotherapy). Resistance to these drugs can emerge due to a variety of reasons including genetic or epigenetic changes which alter the binding site of the drug target, cellular metabolism or export mechanisms. Obtaining a better understanding of the evolution of resistant populations during therapy may enable the design of more effective therapeutic regimens which prevent or delay progression of disease due to resistance.

Evolution of resistance to targeted anti-cancer therapies during continuous and pulsed administration strategies. One important question in clinical cancer research is the identification of optimum therapeutic administration strategies so that the risk of resistance is minimized. In this work we investigated optimal drug dosing schedules to prevent, or at least delay, the emergence of resistance. We designed a stochastic mathematical model describing the evolutionary dynamics of drug resistance emerging in tumor cell population during pulsed and metronomic anti-cancer therapy. We considered drug resistance emerging due to a single (epi)genetic alteration and calculate the probability of resistance arising during specific dosing strategies. We then optimize treatment protocols such that the risk of resistance is minimal while considering drug toxicity and side effects as constraints. This methodology can be used to identify optimum drug administration schedules to avoid resistance conferred by one (epi)genetic alteration for general cancers and treatment types. Link (w/F. Michor)

Evolution of resistance to anti-cancer therapy during general dosing schedules. In previous work we studied the dynamics of resistance during pulsed and continuous administration strategies of targeted drugs. We built upon this framework to study the effect of time-varying dosing schedules and pharmacokinetic effects on the evolutionary dynamics of resistance. The populations of sensitive and resistant cells are modeled as multi-type non-homogeneous birth-death processes in which the drug concentration affects the birth and death rates of both the sensitive and resistant cell populations. Thus, the birth and death rate of each cell population is modeled as a general time-dependent function. This extension allows us to consider the effects of more generalized treatment strategies as well as detailed pharmacokinetic phenomena such as drug elimination and accumulation over multiple doses. We develop estimates for the probability of developing resistance and moments of the resistant cell population size. With these estimates, we can optimize treatment schedules over a subspace of tolerated schedules to minimize the risk of disease progression due to resistance as well as locate ideal schedules for controlling the population size of resistant clones in situations where resistance is inevitable. PDF (w/F. Michor)

Diversity in pre-existing resistance to BCR-ABL inhibitors in Chronic Myeloid Leukemia. Despite the successes of imatinib in treating CML, drug resistance leads to reduced efficacy and progression of disease. Several point mutations in the BCR-ABL kinase domain have been identified that confer resistance to imatinib and second generation BCR-ABL inhibitors dasatinib and nilotinib. The existence and characteristics of resistant cells at the start of therapy are difficult to ascertain, since these clones may be present at frequencies below the detection limit of genomic profiling methods. In this paper, we combine a theoretical modeling approach with recent experimental data to investigate the likelihood, composition, and diversity of pre-existing resistance as well the impact of these factors on the response to targeted therapies. We find that in most patients, there will be at most one resistant cell type at detection; the likelihood of more than one resistant type present at detection is small unless the disease is detected late. We also find that the number of patients that harbor the T315I mutation at detection is similar to the number of patients harboring any other resistant mutation; however, the T315I population is more likely to establish a large-sized clone of resistant cells at detection than other types. Delayed detection increases the risk of pre-existing resistance, and also increases the diversity of the resistant population at diagnosis. Lastly, the overall benefit of combination therapy over monotherapy with imatinib or dasatinib is significantly more pronounced for patients whose disease is diagnosed late. (w/ K. Leder, B. Skaggs, M. Gorre, F. Michor)

Numerical methods for stochastic partial differential equations

Cellular signaling networks, often modeled as high-dimensional PDE or ODE systems, describe the complex signaling interaction structures governing basic cellular activities which form the basis of development, tissue repair, and many other vital processes. Stochasticity in these often hypersensitive networks can arise from many sources, including but not limited to genetic heterogeneity. I am interested in developing tools to study how heterogeneity and variability at the individual cellular level can infuence the outcome of these biochemical signaling networks. More broadly, I am interested in the development of stochastic spectral methods which can be used to quantify the effect of random sources in a PDE or ODE system as it propagates in time via estimation of moments of functionals of the system solutions. These techniques can be applied to study the propagation of noise in many types of physical and biological systems.

The Multi-Element Probabilistic Collocation Method (ME-PCM): Error Analysis and Applications. Stochastic spectral methods are numerical techniques for approximating solutions to partial differential equations with random parameters. In this work, we present and examine the multi-element probabilistic collocation method (ME-PCM), which is a generalized form of the probabilistic collocation method. In the ME-PCM, the parametric space is discretized and a collocation/cubature grid is prescribed on each element. Both full and sparse tensor product grids based on Gauss and Clenshaw- Curtis quadrature rules are considered. We prove analytically and observe in numerical tests that as the parameter space mesh is refined, the convergence rate of the solution depends on the quadrature rule of each element only through its degree of exactness. In addition, the L2 error of the tensor product interpolant is exam- ined and an adaptivity algorithm is provided. Numerical examples demonstrating adaptive ME-PCM are shown, including low-regularity problems and long-time in- tegration. We test the ME-PCM on two-dimensional Navier Stokes examples and a stochastic diffusion problem with various random input distributions and up to 50 dimensions. While the convergence rate of ME-PCM deteriorates in 50 dimensions, the error in the mean and variance is two orders of magnitude lower than the er- ror obtained with the Monte Carlo method using only a small number of samples (e.g., 100). The computational cost of ME-PCM is found to be favorable when com- pared to the cost of other methods including stochastic Galerkin, Monte Carlo and quasi-random sequence methods. PDF (w/X. Wan, G. Karniadakis)

Multi-Element Probabilistic Collocation for Sensitivity Analysis in Cellular Signaling Networks. In this work we formulate the Multi-Element Probabilistic Collocation Method (ME-PCM) as a tool for sensitivity analysis of differential equation models as applied to cellular signaling networks. This method utilizes a simple, efficient sampling algorithm to quantify local sensitivi- ties throughout the parameter space. We apply the ME-PCM to a previously published ordinary differential equation model of the apoptosis signaling network. We first verify agreement with the previously identified regions of sensitivity and then go on to analyze this region in greater detail with the ME-PCM. We demonstrate the generality of the ME-PCM by studying sensitiv- ity of the system using a variety of biologically relevant markers in the system such as variation in one (or many) chemical species as a function of time, and total exposure of a single chemical species. PDF (w/S. Sindi, G. Karniadakis)

Uncertainty quantification for PDEs with high dimensional random inputs. Recently there has been much interest in developing stochastic spectral methods to deal with high-dimensional inputs. I have been working on combining multi-element polynomial chaos with an analysis of variance (ANOVA) functional decomposition to enhance the convergence rate of polynomial chaos in high dimensions and in problems with low stochastic regularity. We have investigated the dependence of the convergence of this method on the ANOVA decomposition parameters and concluded that for problems with prohibitively high nominal dimension and possibly low stochastic regularity, this method is a useful tool in solving problems with random inputs with low effective dimension or a decay in interaction weights. PDF (w/ G. Karniadakis)

Algorithm refinement for the stochastic Burgers equation. In this work, we developed an algorithm refinement (AR) scheme for an excluded random walk model whose mean field behavior is given by the viscous Burgers equation. AR hybrids use the adaptive mesh refinement framework to model a system using a molecular algorithm where desired while allowing a computationally faster continuum representation to be used in the remainder of the domain. The focus in this work is the role of fluctuations on the dynamics. In particular, we demonstrate that it is necessary to include a stochastic forcing term in Burgers equation to accurately capture the correct behavior of the system. The conclusion we draw from this study is that the fidelity of multiscale methods that couple disparate algorithms depends on the consistent modeling of fluctuations in each algorithm and on a coupling, such as algorithm refinement, that preserves this consistency. PDF (w/ J. Bell, A. Garcia)

Stochastic simulation of riser-sections with uncertain measured pressure loads and/or uncertain material properties We investigate three-dimensional problems in solid mechanics with stochastic loading or material properties. To solve these problems, we use a spectral expansion of the solution and random inputs based on Askey-type orthogonal polynomials in terms of independent, identically distributed (i.i.d) random variables. A Galerkin procedure using these types of expansions, the generalized Polynomial Chaos (gPC) method, is employed to solve linear elasticity problems. An analagous spectral collocation formulation is used to study problems in nonlinear elasticity. These methods both cast the stochastic problem as a coupled or decoupled high-dimensional system of deterministic PDEs, which is then solved numerically using a deterministic p-finite element solver. We present algorithms for solving certain coupled systems arising from the stochastic Galerkin projection without modifying the original deterministic solver. Three-dimensional riser-sections undergoing elastic deformations due to random pressure loads are considered. We also model a riser-section with stochastic Youngs modulus undergoing deterministic loads. It is demonstrated that the gPC method provides accurate and efficient results at a speed-up factor of two and three orders of magnitude compared to traditional Monte-Carlo simulations. For nonlinear problems, the stochastic collocation method is also shown to be much faster than Monte-Carlo simulation, while still rivaling this method in simplicity of implementation. PDF (w/ Z. Yosibash, G. Karniadakis)