A Multi-stage Representation of Cell Proliferation as a Markov Process

A comprehensive review of computational cell cycle models in guiding cancer treatment strategies

This article reviews the current knowledge and recent advancements in computational modeling of the cell cycle. It offers a comparative analysis of various modeling paradigms, highlighting their unique strengths, limitations, and applications. Specifically, the article compares deterministic and stochastic models, single-cell versus population models, and mechanistic versus abstract models. This detailed analysis helps determine the most suitable modeling framework for various research needs. Additionally, the discussion extends to the utilization of these computational models to illuminate cell cycle dynamics, with a particular focus on cell cycle viability, crosstalk with signaling pathways, tumor microenvironment, DNA replication, and repair mechanisms, underscoring their critical roles in tumor progression and the optimization of cancer therapies. By applying these models to crucial aspects of cancer therapy planning for better outcomes, including drug efficacy quantification, drug discovery, drug resistance analysis, and dose optimization, the review highlights the significant potential of computational insights in enhancing the precision and effectiveness of cancer treatments. This emphasis on the intricate relationship between computational modeling and therapeutic strategy development underscores the pivotal role of advanced modeling techniques in navigating the complexities of cell cycle dynamics and their implications for cancer therapy.

Quantifying cell cycle regulation by tissue crowding

The spatiotemporal coordination and regulation of cell proliferation is fundamental in many aspects of development and tissue maintenance. Cells have the ability to adapt their division rates in response to mechanical constraints, yet we do not fully understand how cell proliferation regulation impacts cell migration phenomena. Here, we present a minimal continuum model of cell migration with cell cycle dynamics, which includes density-dependent effects and hence can account for cell proliferation regulation. By combining minimal mathematical modeling, Bayesian inference, and recent experimental data, we quantify the impact of tissue crowding across different cell cycle stages in epithelial tissue expansion experiments. Our model suggests that cells sense local density and adapt cell cycle progression in response, during G1 and the combined S/G2/M phases, providing an explicit relationship between each cell-cycle-stage duration and local tissue density, which is consistent with several experimental observations. Finally, we compare our mathematical model's predictions to different experiments studying cell cycle regulation and present a quantitative analysis on the impact of density-dependent regulation on cell migration patterns. Our work presents a systematic approach for investigating and analyzing cell cycle data, providing mechanistic insights into how individual cells regulate proliferation, based on population-based experimental measurements.

Stochastic model of Alzheimer's disease progression using two-state Markov chains

In 2016, Hao and Friedman developed a deterministic model of Alzheimer's disease progression using a system of partial differential equations. This model describes the general behavior of the disease, however, it does not incorporate the molecular and cellular stochasticity intrinsic to the underlying disease processes. Here we extend the Hao and Friedman model by modeling each event in disease progression as a stochastic Markov process. This model identifies stochasticity in disease progression, as well as changes to the mean dynamics of key agents. We find that the pace of neuron death increases whereas the production of the two key measures of progression, Tau and Amyloid beta proteins, decelerates when stochasticity is incorporated into the model. These results suggest that the non-constant reactions and time-steps have a significant effect on the overall progression of the disease.

Mathematical Modeling Identifies Optimum Palbociclib-fulvestrant Dose Administration Schedules for the Treatment of Patients with Estrogen Receptor-positive Breast Cancer

Cyclin-dependent kinases 4/6 (CDK4/6) inhibitors such as palbociclib are approved for the treatment of metastatic estrogen receptor-positive (ER+) breast cancer in combination with endocrine therapies and significantly improve outcomes in patients with this disease. However, given the large number of possible pairwise drug combinations and administration schedules, it remains unclear which clinical strategy would lead to best survival. Here, we developed a computational, cell cycle-explicit model to characterize the pharmacodynamic response to palbociclib-fulvestrant combination therapy. This pharmacodynamic model was parameterized, in a Bayesian statistical inference approach, using in vitro data from cells with wild-type estrogen receptor (WT-ER) and cells expressing the activating missense ER mutation, Y537S, which confers resistance to fulvestrant. We then incorporated pharmacokinetic models derived from clinical data into our computational modeling platform. To systematically compare dose administration schedules, we performed in silico clinical trials based on integrating our pharmacodynamic and pharmacokinetic models as well as considering clinical toxicity constraints. We found that continuous dosing of palbociclib is more effective for lowering overall tumor burden than the standard, pulsed-dose palbociclib treatment. Importantly, our mathematical modeling and statistical analysis platform provides a rational method for comparing treatment strategies in search of optimal combination dosing strategies of other cell-cycle inhibitors in ER+ breast cancer.

SIGNIFICANCE

We created a computational modeling platform to predict the effects of fulvestrant/palbocilib treatment on WT-ER and Y537S-mutant breast cancer cells, and found that continuous treatment schedules are more effective than the standard, pulsed-dose palbociclib treatment schedule.

Distribution modeling quantifies collective TH cell decision circuits in chronic inflammation

Immune responses are tightly regulated by a diverse set of interacting immune cell populations. Alongside decision-making processes such as differentiation into specific effector cell types, immune cells initiate proliferation at the beginning of an inflammation, forming two layers of complexity. Here, we developed a general mathematical framework for the data-driven analysis of collective immune cell dynamics. We identified qualitative and quantitative properties of generic network motifs, and we specified differentiation dynamics by analysis of kinetic transcriptome data. Furthermore, we derived a specific, data-driven mathematical model for T helper 1 versus T follicular helper cell-fate decision dynamics in acute and chronic lymphocytic choriomeningitis virus infections in mice. The model recapitulates important dynamical properties without model fitting and solely by using measured response-time distributions. Model simulations predict different windows of opportunity for perturbation in acute and chronic infection scenarios, with potential implications for optimization of targeted immunotherapy.

Fluctuations in tissue growth portray homeostasis as a critical state and long-time non-Markovian cell proliferation as Markovian

Tissue growth is an emerging phenomenon that results from the cell-level interplay between proliferation and apoptosis, which is crucial during embryonic development, tissue regeneration, as well as in pathological conditions such as cancer. In this theoretical article, we address the problem of stochasticity in tissue growth by first considering a minimal Markovian model of tissue size, quantified as the number of cells in a simulated tissue, subjected to both proliferation and apoptosis. We find two dynamic phases, growth and decay, separated by a critical state representing a homeostatic tissue. Since the main limitation of the Markovian model is its neglect of the cell cycle, we incorporated a refractory period that temporarily prevents proliferation immediately following cell division, as a minimal proxy for the cell cycle, and studied the model in the growth phase. Importantly, we obtained from this last model an effective Markovian rate, which accurately describes general trends of tissue size. This study shows that the dynamics of tissue growth can be theoretically conceptualized as a Markovian process where homeostasis is a critical state flanked by decay and growth phases. Notably, in the growing non-Markovian model, a Markovian-like growth process emerges at large time scales.

Solving the stochastic dynamics of population growth

Population growth is a fundamental process in ecology and evolution. The population size dynamics during growth are often described by deterministic equations derived from kinetic models. Here, we simulate several population growth models and compare the size averaged over many stochastic realizations with the deterministic predictions. We show that these deterministic equations are generically bad predictors of the average stochastic population dynamics. Specifically, deterministic predictions overestimate the simulated population sizes, especially those of populations starting with a small number of individuals. Describing population growth as a stochastic birth process, we prove that the discrepancy between deterministic predictions and simulated data is due to unclosed-moment dynamics. In other words, the deterministic approach does not consider the variability of birth times, which is particularly important with small population sizes. We show that some moment-closure approximations describe the growth dynamics better than the deterministic prediction. However, they do not reduce the error satisfactorily and only apply to some population growth models. We explicitly solve the stochastic growth dynamics, and our solution applies to any population growth model. We show that our solution exactly quantifies the dynamics of a community composed of different strains and correctly predicts the fixation probability of a strain in a serial dilution experiment. Our work sets the foundations for a more faithful modeling of community and population dynamics. It will allow the development of new tools for a more accurate analysis of experimental and empirical results, including the inference of important growth parameters.

Stochastic model of Alzheimer's Disease progression using two-state Markov chains

In 2016, Hao and Friedman developed a deterministic model of Alzheimer's disease progression using a system of partial differential equations. This model describes the general behavior of the disease, however, it does not incorporate the molecular and cellular stochasticity intrinsic to the underlying disease processes. Here we extend the Hao and Friedman model by modeling each event in disease progression as a stochastic Markov process. This model identifies stochasticity in disease progression, as well as changes to the mean dynamics of key agents. We find that the pace of neuron death increases whereas the production of the two key measures of progression, Tau and Amyloid beta proteins, decelerates when stochasticity is incorporated into the model. These results suggest that the non-constant reactions and time-steps have a significant effect on the overall progression of the disease.

Analysis and modeling of cancer drug responses using cell cycle phase-specific rate effects

Identifying effective therapeutic treatment strategies is a major challenge to improving outcomes for patients with breast cancer. To gain a comprehensive understanding of how clinically relevant anti-cancer agents modulate cell cycle progression, here we use genetically engineered breast cancer cell lines to track drug-induced changes in cell number and cell cycle phase to reveal drug-specific cell cycle effects that vary across time. We use a linear chain trick (LCT) computational model, which faithfully captures drug-induced dynamic responses, correctly infers drug effects, and reproduces influences on specific cell cycle phases. We use the LCT model to predict the effects of unseen drug combinations and confirm these in independent validation experiments. Our integrated experimental and modeling approach opens avenues to assess drug responses, predict effective drug combinations, and identify optimal drug sequencing strategies.

On predicting heterogeneity in nanoparticle dosage

Nanoparticles are increasingly employed as a vehicle for the targeted delivery of therapeutics to specific cell types. However, much remains to be discovered about the fundamental biology that dictates the interactions between nanoparticles and cells. Accordingly, few nanoparticle-based targeted therapeutics have succeeded in clinical trials. One element that hinders our understanding of nanoparticle-cell interactions is the presence of heterogeneity in nanoparticle dosage data obtained from standard experiments. It is difficult to distinguish between heterogeneity that arises from stochasticity in nanoparticle-cell interactions, and that which arises from heterogeneity in the cell population. Mathematical investigations have revealed that both sources of heterogeneity contribute meaningfully to the heterogeneity in nanoparticle dosage. However, these investigations have relied on simplified models of nanoparticle internalisation. Here we present a stochastic mathematical model of nanoparticle internalisation that incorporates a suite of relevant biological phenomena such as multistage internalisation, cell division, asymmetric nanoparticle inheritance and nanoparticle saturation. Critically, our model provides information about nanoparticle dosage at an individual cell level. We perform model simulations to examine the influence of specific biological phenomena on the heterogeneity in nanoparticle dosage in the absence of heterogeneity in the cell population. Under certain modelling assumptions, we derive analytic approximations of the nanoparticle dosage distribution. We demonstrate that the analytic approximations are accurate, and show that nanoparticle dosage can be described by a Poisson mixture distribution with rate parameters that are a function of Beta-distributed random variables. We discuss the implications of the analytic results with respect to parameter estimation and model identifiability from standard experimental data. Finally, we highlight extensions and directions for future research.

Concentration fluctuations in growing and dividing cells: Insights into the emergence of concentration homeostasis

Intracellular reaction rates depend on concentrations and hence their levels are often regulated. However classical models of stochastic gene expression lack a cell size description and cannot be used to predict noise in concentrations. Here, we construct a model of gene product dynamics that includes a description of cell growth, cell division, size-dependent gene expression, gene dosage compensation, and size control mechanisms that can vary with the cell cycle phase. We obtain expressions for the approximate distributions and power spectra of concentration fluctuations which lead to insight into the emergence of concentration homeostasis. We find that (i) the conditions necessary to suppress cell division-induced concentration oscillations are difficult to achieve; (ii) mRNA concentration and number distributions can have different number of modes; (iii) two-layer size control strategies such as sizer-timer or adder-timer are ideal because they maintain constant mean concentrations whilst minimising concentration noise; (iv) accurate concentration homeostasis requires a fine tuning of dosage compensation, replication timing, and size-dependent gene expression; (v) deviations from perfect concentration homeostasis show up as deviations of the concentration distribution from a gamma distribution. Some of these predictions are confirmed using data for E. coli, fission yeast, and budding yeast.

Counting generations in birth and death processes with competing Erlang and exponential waiting times

Lymphocyte populations, stimulated in vitro or in vivo, grow as cells divide. Stochastic models are appropriate because some cells undergo multiple rounds of division, some die, and others of the same type in the same conditions do not divide at all. If individual cells behave independently, then each cell can be imagined as sampling from a probability density of times to division and death. The exponential density is the most mathematically and computationally convenient choice. It has the advantage of satisfying the memoryless property, consistent with a Markov process, but it overestimates the probability of short division times. With the aim of preserving the advantages of a Markovian framework while improving the representation of experimentally-observed division times, we consider a multi-stage model of cellular division and death. We use Erlang-distributed (or, more generally, phase-type distributed) times to division, and exponentially distributed times to death. We classify cells into generations, using the rule that the daughters of cells in generation n are in generation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n+1$$\end{document}n+1. In some circumstances, our representation is equivalent to established models of lymphocyte dynamics. We find the growth rate of the cell population by calculating the proportions of cells by stage and generation. The exponent describing the late-time cell population growth, and the criterion for extinction of the population, differs from what would be expected if N steps with rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document}λ were equivalent to a single step of rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda /N$$\end{document}λ/N. We link with a published experimental dataset, where cell counts were reported after T cells were transferred to lymphopenic mice, using Approximate Bayesian Computation. In the comparison, the death rate is assumed to be proportional to the generation and the Erlang time to division for generation 0 is allowed to differ from that of subsequent generations. The multi-stage representation is preferred to a simple exponential in posterior distributions, and the mean time to first division is estimated to be longer than the mean time to subsequent divisions.

Equivalence framework for an age-structured multistage representation of the cell cycle

We develop theoretical equivalences between stochastic and deterministic models for populations of individual cells stratified by age. Specifically, we develop a hierarchical system of equations describing the full dynamics of an age-structured multistage Markov process for approximating cell cycle time distributions. We further demonstrate that the resulting mean behavior is equivalent, over large timescales, to the classical McKendrick-von Foerster integropartial differential equation. We conclude by extending this framework to a spatial context, facilitating the modeling of traveling wave phenomena and cell-mediated pattern formation. More generally, this methodology may be extended to myriad reaction-diffusion processes for which the age of individuals is relevant to the dynamics.

Cyton2: A Model of Immune Cell Population Dynamics That Includes Familial Instructional Inheritance

Lymphocytes are the central actors in adaptive immune responses. When challenged with antigen, a small number of B and T cells have a cognate receptor capable of recognising and responding to the insult. These cells proliferate, building an exponentially growing, differentiating clone army to fight off the threat, before ceasing to divide and dying over a period of weeks, leaving in their wake memory cells that are primed to rapidly respond to any repeated infection. Due to the non-linearity of lymphocyte population dynamics, mathematical models are needed to interrogate data from experimental studies. Due to lack of evidence to the contrary and appealing to arguments based on Occam’s Razor, in these models newly born progeny are typically assumed to behave independently of their predecessors. Recent experimental studies, however, challenge that assumption, making clear that there is substantial inheritance of timed fate changes from each cell by its offspring, calling for a revision to the existing mathematical modelling paradigms used for information extraction. By assessing long-term live-cell imaging of stimulated murine B and T cells in vitro, we distilled the key phenomena of these within-family inheritances and used them to develop a new mathematical model, Cyton2, that encapsulates them. We establish the model’s consistency with these newly observed fine-grained features. Two natural concerns for any model that includes familial correlations would be that it is overparameterised or computationally inefficient in data fitting, but neither is the case for Cyton2. We demonstrate Cyton2’s utility by challenging it with high-throughput flow cytometry data, which confirms the robustness of its parameter estimation as well as its ability to extract biological meaning from complex mixed stimulation experiments. Cyton2, therefore, offers an alternate mathematical model, one that is, more aligned to experimental observation, for drawing inferences on lymphocyte population dynamics.

αβ/γδ T cell lineage outcome is regulated by intrathymic cell localization and environmental signals

αβ and γδ T cells are two distinct sublineages that develop in the vertebrate thymus. Thus far, their differentiation from a common progenitor is mostly understood to be regulated by intrinsic mechanisms. However, the proportion of αβ/γδ T cells varies in different vertebrate taxa. How this process is regulated in species that tend to produce a high frequency of γδ T cells is unstudied. Using an in vivo teleost model, the medaka, we report that progenitors first enter a thymic niche where their development into γδ T cells is favored. Translocation from this niche, mediated by chemokine receptor Ccr9b, is a prerequisite for their differentiation into αβ T cells. On the other hand, the thymic niche also generates opposing gradients of the cytokine interleukin-7 and chemokine Ccl25a, and, together, they influence the lineage outcome. We propose a previously unknown mechanism that determines the proportion of αβ/γδ lineages within species.

Building mean field ODE models using the generalized linear chain trick & Markov chain theory

The well known linear chain trick (LCT) allows modellers to derive mean field ODEs that assume gamma (Erlang) distributed passage times, by transitioning individuals sequentially through a chain of sub-states. The time spent in these sub-states is the sum of k exponentially distributed random variables, and is thus gamma distributed. The generalized linear chain trick (GLCT) extends this technique to the broader phase-type family of distributions, which includes exponential, Erlang, hypoexponential, and Coxian distributions. Phase-type distributions are the family of matrix exponential distributions on [0,∞) that represent the absorption time distributions for finite-state, continuous time Markov chains (CTMCs). Here we review CTMCs and phase-type distributions, then illustrate how to use the GLCT to efficiently build ODE models from underlying stochastic model assumptions. We introduce two novel model families by using the GLCT to generalize the Rosenzweig-MacArthur predator-prey model, and the SEIR model. We illustrate the kinds of complexity that can be captured by such models through multiple examples. We also show the benefits of using a GLCT-based model formulation to speed up the computation of numerical solutions to such models. These results highlight the intuitive nature, and utility, of using the GLCT to derive ODE models from first principles.

A mathematical model of viral oncology as an immuno-oncology instigator

We develop and analyse a mathematical model of tumour-immune interaction that explicitly incorporates heterogeneity in tumour cell cycle duration by using a distributed delay differential equation. We derive a necessary and sufficient condition for local stability of the cancer-free equilibrium in which the amount of tumour-immune interaction completely characterizes disease progression. Consistent with the immunoediting hypothesis, we show that decreasing tumour-immune interaction leads to tumour expansion. Finally, by simulating the mathematical model, we show that the strength of tumour-immune interaction determines the long-term success or failure of viral therapy.

Optimal transition paths of phenotypic switching in a non-Markovian self-regulation gene expression

Gene expression is a complex biochemical process involving multiple reaction steps, creating molecular memory because the probability of waiting time between consecutive reaction steps no longer follows exponential distributions. What effect the molecular memory has on metastable states in gene expression remains not fully understood. Here, we study transition paths of switching between bistable states for a non-Markovian model of gene expression equipped with a self-regulation. Employing the large deviation theory for this model, we analyze the optimal transition paths of switching between bistable states in gene expression, interestingly finding that dynamic behaviors in gene expression along the optimal transition paths significantly depend on the molecular memory. Moreover, we discover that the molecular memory can prolong the time of switching between bistable states in gene expression along the optimal transition paths. Our results imply that the molecular memory may be an unneglectable factor to affect switching between metastable states in gene expression.

Cell size distribution of lineage data: analytic results and parameter inference

Recent advances in single-cell technologies have enabled time-resolved measurements of the cell size over several cell cycles. These data encode information on how cells correct size aberrations so that they do not grow abnormally large or small. Here, we formulate a piecewise deterministic Markov model describing the evolution of the cell size over many generations, for all three cell size homeostasis strategies (timer, sizer, and adder). The model is solved to obtain an analytical expression for the non-Gaussian cell size distribution in a cell lineage; the theory is used to understand how the shape of the distribution is influenced by the parameters controlling the dynamics of the cell cycle and by the choice of cell tracking protocol. The theoretical cell size distribution is found to provide an excellent match to the experimental cell size distribution of E. coli lineage data collected under various growth conditions.

A novel mathematical model of heterogeneous cell proliferation

We present a novel mathematical model of heterogeneous cell proliferation where the total population consists of a subpopulation of slow-proliferating cells and a subpopulation of fast-proliferating cells. The model incorporates two cellular processes, asymmetric cell division and induced switching between proliferative states, which are important determinants for the heterogeneity of a cell population. As motivation for our model we provide experimental data that illustrate the induced-switching process. Our model consists of a system of two coupled delay differential equations with distributed time delays and the cell densities as functions of time. The distributed delays are bounded and allow for the choice of delay kernel. We analyse the model and prove the nonnegativity and boundedness of solutions, the existence and uniqueness of solutions, and the local stability characteristics of the equilibrium points. We find that the parameters for induced switching are bifurcation parameters and therefore determine the long-term behaviour of the model. Numerical simulations illustrate and support the theoretical findings, and demonstrate the primary importance of transient dynamics for understanding the evolution of many experimental cell populations.

Modeling low-dose radiation-induced acute myeloid leukemia in male CBA/H mice

The effect of low-dose ionizing radiation exposure on leukemia incidence remains poorly understood. Possible dose-response curves for various forms of leukemia are largely based on cohorts of atomic bomb survivors. Animal studies can contribute to an improved understanding of radiation-induced acute myeloid leukemia (rAML) in humans. In male CBA/H mice, incidence of rAML can be described by a two-hit model involving a radiation-induced deletion with Sfpi1 gene copy loss and a point mutation in the remaining Sfpi1 allele. In the present study (historical) mouse data were used and these processes were translated into a mathematical model to study photon-induced low-dose AML incidence in male CBA/H mice following acute exposure. Numerical model solutions for low-dose rAML incidence and diagnosis times could respectively be approximated with a model linear-quadratic in radiation dose and a normal cumulative distribution function. Interestingly, the low-dose incidence was found to be proportional to the modeled number of cells carrying the Sfpi1 deletion present per mouse following exposure. After making only model-derived high-dose rAML estimates available to extrapolate from, the linear-quadratic model could be used to approximate low-dose rAML incidence calculated with our mouse model. The accuracy in estimating low-dose rAML incidence when extrapolating from a linear model using a low-dose effectiveness factor was found to depend on whether a data transformation was used in the curve fitting procedure.

Synchronized oscillations in growing cell populations are explained by demographic noise

Understanding synchrony in growing populations is important for applications as diverse as epidemiology and cancer treatment. Recent experiments employing fluorescent reporters in melanoma cell lines have uncovered growing subpopulations exhibiting sustained oscillations, with nearby cells appearing to synchronize their cycles. In this study, we demonstrate that the behavior observed is consistent with long-lasting transient phenomenon initiated and amplified by the finite-sample effects and demographic noise. We present a novel mathematical analysis of a multistage model of cell growth, which accurately reproduces the synchronized oscillations. As part of the analysis, we elucidate the transient and asymptotic phases of the dynamics and derive an analytical formula to quantify the effect of demographic noise in the appearance of the oscillations. The implications of these findings are broad, such as providing insight into experimental protocols that are used to study the growth of asynchronous populations and, in particular, those investigations relating to anticancer drug discovery.

Heterogeneity and 'memory' in stem cell populations

Modern single cell experiments have revealed unexpected heterogeneity in apparently functionally 'pure' cell populations. However, we are still lacking a conceptual framework to understand this heterogeneity. Here, we propose that cellular memories-changes in the molecular status of a cell in response to a stimulus, that modify the ability of the cell to respond to future stimuli-are an essential ingredient in any such theory. We illustrate this idea by considering a simple age-structured model of stem cell proliferation that takes account of mitotic memories. Using this model we argue that asynchronous mitosis generates heterogeneity that is central to stem cell population function. This model naturally explains why stem cell numbers increase through life, yet regenerative potency simultaneously declines.

Asymmetry of nanoparticle inheritance upon cell division: Effect on the coefficient of variation

Several previous studies have shown that when a cell that has taken up nanoparticles divides, the nanoparticles are inherited by the two daughter cells in an asymmetrical fashion, with one daughter cell receiving more nanoparticles than the other. This interesting observation is typically demonstrated either indirectly using mathematical modelling of high-throughput experimental data or more directly by imaging individual cells as they divide. Here we suggest that measurements of the coefficient of variation (standard deviation over mean) of the number of nanoparticles per cell over the cell population is another means of assessing the degree of asymmetry. Using simulations of an evolving cell population, we show that the coefficient of variation is sensitive to the degree of asymmetry and note its characteristic evolution in time. As the coefficient of variation is readily measurable using high-throughput techniques, this should allow a more rapid experimental assessment of the degree of asymmetry.

An in silico model of LINE-1-mediated neoplastic evolution

MOTIVATION

Recent research has uncovered roles for transposable elements (TEs) in multiple evolutionary processes, ranging from somatic evolution in cancer to putatively adaptive germline evolution across species. Most models of TE population dynamics, however, have not incorporated actual genome sequence data. The effect of site integration preferences of specific TEs on evolutionary outcomes and the effects of different selection regimes on TE dynamics in a specific genome are unknown. We present a stochastic model of LINE-1 (L1) transposition in human cancer. This system was chosen because the transposition of L1 elements is well understood, the population dynamics of cancer tumors has been modeled extensively, and the role of L1 elements in cancer progression has garnered interest in recent years.

RESULTS

Our model predicts that L1 retrotransposition (RT) can play either advantageous or deleterious roles in tumor progression, depending on the initial lesion size, L1 insertion rate and tumor driver genes. Small changes in the RT rate or set of driver tumor-suppressor genes (TSGs) were observed to alter the dynamics of tumorigenesis. We found high variation in the density of L1 target sites across human protein-coding genes. We also present an analysis, across three cancer types, of the frequency of homozygous TSG disruption in wild-type hosts compared to those with an inherited driver allele.

AVAILABILITY AND IMPLEMENTATION

Source code is available at https://github.com/atallah-lab/neoplastic-evolution.

CONTACT

jlebien@uno.edu.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Effects of cell cycle variability on lineage and population measurements of messenger RNA abundance

Many models of gene expression do not explicitly incorporate a cell cycle description. Here, we derive a theory describing how messenger RNA (mRNA) fluctuations for constitutive and bursty gene expression are influenced by stochasticity in the duration of the cell cycle and the timing of DNA replication. Analytical expressions for the moments show that omitting cell cycle duration introduces an error in the predicted mean number of mRNAs that is a monotonically decreasing function of η, which is proportional to the ratio of the mean cell cycle duration and the mRNA lifetime. By contrast, the error in the variance of the mRNA distribution is highest for intermediate values of η consistent with genome-wide measurements in many organisms. Using eukaryotic cell data, we estimate the errors in the mean and variance to be at most 3% and 25%, respectively. Furthermore, we derive an accurate negative binomial mixture approximation to the mRNA distribution. This indicates that stochasticity in the cell cycle can introduce fluctuations in mRNA numbers that are similar to the effect of bursty transcription. Finally, we show that for real experimental data, disregarding cell cycle stochasticity can introduce errors in the inference of transcription rates larger than 10%.

Isolating the sources of heterogeneity in nano-engineered particle-cell interactions

Nano-engineered particles have the potential to enhance therapeutic success and reduce toxicity-based treatment side effects via the targeted delivery of drugs to cells. This delivery relies on complex interactions between numerous biological, chemical and physical processes. The intertwined nature of these processes has thus far hindered attempts to understand their individual impact. Variation in experimental data, such as the number of particles inside each cell, further inhibits understanding. Here, we present a mathematical framework that is capable of examining the impact of individual processes during particle delivery. We demonstrate that variation in experimental particle uptake data can be explained by three factors: random particle motion; variation in particle-cell interactions; and variation in the maximum particle uptake per cell. Without all three factors, the experimental data cannot be explained. This work provides insight into biological mechanisms that cause heterogeneous responses to treatment, and enables precise identification of treatment-resistant cell subpopulations.

Exact solution of stochastic gene expression models with bursting, cell cycle and replication dynamics

The bulk of stochastic gene expression models in the literature do not have an explicit description of the age of a cell within a generation and hence they cannot capture events such as cell division and DNA replication. Instead, many models incorporate the cell cycle implicitly by assuming that dilution due to cell division can be described by an effective decay reaction with first-order kinetics. If it is further assumed that protein production occurs in bursts, then the stationary protein distribution is a negative binomial. Here we seek to understand how accurate these implicit models are when compared with more detailed models of stochastic gene expression. We derive the exact stationary solution of the chemical master equation describing bursty protein dynamics, binomial partitioning at mitosis, age-dependent transcription dynamics including replication, and random interdivision times sampled from Erlang or more general distributions; the solution is different for single lineage and population snapshot settings. We show that protein distributions are well approximated by the solution of implicit models (a negative binomial) when the mean number of mRNAs produced per cycle is low and the cell cycle length variability is large. When these conditions are not met, the distributions are either almost bimodal or else display very flat regions near the mode and cannot be described by implicit models. We also show that for genes with low transcription rates, the size of protein noise has a strong dependence on the replication time, it is almost independent of cell cycle variability for lineage measurements, and increases with cell cycle variability for population snapshot measurements. In contrast for large transcription rates, the size of protein noise is independent of replication time and increases with cell cycle variability for both lineage and population measurements.

Analytical distributions for detailed models of stochastic gene expression in eukaryotic cells

The stochasticity of gene expression presents significant challenges to the modeling of genetic networks. A two-state model describing promoter switching, transcription, and messenger RNA (mRNA) decay is the standard model of stochastic mRNA dynamics in eukaryotic cells. Here, we extend this model to include mRNA maturation, cell division, gene replication, dosage compensation, and growth-dependent transcription. We derive expressions for the time-dependent distributions of nascent mRNA and mature mRNA numbers, provided two assumptions hold: 1) nascent mRNA dynamics are much faster than those of mature mRNA; and 2) gene-inactivation events occur far more frequently than gene-activation events. We confirm that thousands of eukaryotic genes satisfy these assumptions by using data from yeast, mouse, and human cells. We use the expressions to perform a sensitivity analysis of the coefficient of variation of mRNA fluctuations averaged over the cell cycle, for a large number of genes in mouse embryonic stem cells, identifying degradation and gene-activation rates as the most sensitive parameters. Furthermore, it is shown that, despite the model's complexity, the time-dependent distributions predicted by our model are generally well approximated by the negative binomial distribution. Finally, we extend our model to include translation, protein decay, and auto-regulatory feedback, and derive expressions for the approximate time-dependent protein-number distributions, assuming slow protein decay. Our expressions enable us to study how complex biological processes contribute to the fluctuations of gene products in eukaryotic cells, as well as allowing a detailed quantitative comparison with experimental data via maximum-likelihood methods.

Estimating growth patterns and driver effects in tumor evolution from individual samples

Tumors accumulate thousands of mutations, and sequencing them has given rise to methods for finding cancer drivers via mutational recurrence. However, these methods require large cohorts and underperform for low recurrence. Recently, ultra-deep sequencing has enabled accurate measurement of VAFs (variant-allele frequencies) for mutations, allowing the determination of evolutionary trajectories. Here, based solely on the VAF spectrum for an individual sample, we report on a method that identifies drivers and quantifies tumor growth. Drivers introduce perturbations into the spectrum, and our method uses the frequency of hitchhiking mutations preceding a driver to measure this. As validation, we use simulation models and 993 tumors from the Pan-Cancer Analysis of Whole Genomes (PCAWG) Consortium with previously identified drivers. Then we apply our method to an ultra-deep sequenced acute myeloid leukemia (AML) tumor and identify known cancer genes and additional driver candidates. In summary, our framework presents opportunities for personalized driver diagnosis using sequencing data from a single individual.

A stochastic model of adult neurogenesis coupling cell cycle progression and differentiation

Long-term tissue homeostasis requires a precise balance between stem cell self-renewal and the generation of differentiated progeny. Recently, it has been shown that in the adult murine brain, neural stem cells (NSCs) divide mostly symmetrically. This finding suggests that the required balance for tissue homeostasis is accomplished at the population level. However, it remains unclear how this balance is enabled. Furthermore, there is experimental evidence that proneural differentiation factors not only promote differentiation, but also cell cycle progression, suggesting a link between the two processes in NSCs. To study the effect of such a link on NSC dynamics, we developed a stochastic model in which stem cells have an intrinsic probability to progress through cell cycle and to differentiate. Our results show that increasing heterogeneity in differentiation probabilities leads to a decreased probability of long-term tissue homeostasis, and that this effect can be compensated when cell cycle progression and differentiation are positively coupled. Using single-cell RNA-Seq profiling of adult NSCs, we found a positive correlation in the expression levels of cell cycle and differentiation markers. Our findings suggest that a coupling between cell cycle progression and differentiation on the cellular level is part of the process that maintains tissue homeostasis in the adult brain.

Mathematical models incorporating a multi-stage cell cycle replicate normally-hidden inherent synchronization in cell proliferation

We present a suite of experimental data showing that cell proliferation assays, prepared using standard methods thought to produce asynchronous cell populations, persistently exhibit inherent synchronization. Our experiments use fluorescent cell cycle indicators to reveal the normally hidden cell synchronization, by highlighting oscillatory subpopulations within the total cell population. These oscillatory subpopulations would never be observed without these cell cycle indicators. On the other hand, our experimental data show that the total cell population appears to grow exponentially, as in an asynchronous population. We reconcile these seemingly inconsistent observations by employing a multi-stage mathematical model of cell proliferation that can replicate the oscillatory subpopulations. Our study has important implications for understanding and improving experimental reproducibility. In particular, inherent synchronization may affect the experimental reproducibility of studies aiming to investigate cell cycle-dependent mechanisms, including changes in migration and drug response.

Generalizations of the 'Linear Chain Trick': incorporating more flexible dwell time distributions into mean field ODE models

In this paper we generalize the Linear Chain Trick (LCT; aka the Gamma Chain Trick) to help provide modelers more flexibility to incorporate appropriate dwell time assumptions into mean field ODEs, and help clarify connections between individual-level stochastic model assumptions and the structure of corresponding mean field ODEs. The LCT is a technique used to construct mean field ODE models from continuous-time stochastic state transition models where the time an individual spends in a given state (i.e., the dwell time) is Erlang distributed (i.e., gamma distributed with integer shape parameter). Despite the LCT’s widespread use, we lack general theory to facilitate the easy application of this technique, especially for complex models. Modelers must therefore choose between constructing ODE models using heuristics with oversimplified dwell time assumptions, using time consuming derivations from first principles, or to instead use non-ODE models (like integro-differential or delay differential equations) which can be cumbersome to derive and analyze. Here, we provide analytical results that enable modelers to more efficiently construct ODE models using the LCT or related extensions. Specifically, we provide (1) novel LCT extensions for various scenarios found in applications, including conditional dwell time distributions; (2) formulations of these LCT extensions that bypass the need to derive ODEs from integral equations; and (3) a novel Generalized Linear Chain Trick (GLCT) framework that extends the LCT to a much broader set of possible dwell time distribution assumptions, including the flexible phase-type distributions which can approximate distributions on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {R}}^+$$\end{document}R+ and can be fit to data.

Evidence that the human cell cycle is a series of uncoupled, memoryless phases

The cell cycle is canonically described as a series of four consecutive phases: G1, S, G2, and M. In single cells, the duration of each phase varies, but the quantitative laws that govern phase durations are not well understood. Using time-lapse microscopy, we found that each phase duration follows an Erlang distribution and is statistically independent from other phases. We challenged this observation by perturbing phase durations through oncogene activation, inhibition of DNA synthesis, reduced temperature, and DNA damage. Despite large changes in durations in cell populations, phase durations remained uncoupled in individual cells. These results suggested that the independence of phase durations may arise from a large number of molecular factors that each exerts a minor influence on the rate of cell cycle progression. We tested this model by experimentally forcing phase coupling through inhibition of cyclin-dependent kinase 2 (CDK2) or overexpression of cyclin D. Our work provides an explanation for the historical observation that phase durations are both inherited and independent and suggests how cell cycle progression may be altered in disease states.

Neural crest migration with continuous cell states

Models of cranial neural crest cell migration in cell-induced (or self-generated) gradients have included a division of labour into leader and follower migratory states, which undergo chemotaxis and contact guidance, respectively. Despite validated utility of these models through experimental perturbation of migration in the chick embryo and gene expression analysis showing relevant heterogeneity at the single cell level, an often raised question has been whether the discrete cell states are necessary, or if a continuum of cell behaviours offers a functionally equivalent description. Here we argue that this picture is supported by recent single-cell transcriptome data. Motivated by this, we implement two versions of a continuous-state model: (1) signal choice and (2) signal combination. We find that the cell population migrates further than in the discrete-state model and than in experimental observations. We further show that the signal combination model, but not the signal choice model, can be successfully adjusted to experimentally plausible regimes by reducing the chemoattractant consumption parameter. Thus we show an equivalently plausible, experimentally motivated, model of neural crest cell migration.

Proliferation characteristics of cells cultured under periodic versus static conditions

In vitro culture models have become an indispensable tool for assessing a vast variety of biological questions in many scientific fields. However, common in vitro cultures are maintained under static conditions, which do not reflect the in vivo situation and create a non-physiological environment. To assess whether the growth characteristics of cells cultured at pulsed-perfused versus static conditions differ, we observed the growth of differentially cultured cells in vitro by life-cell time-lapse imaging of recombinant HEK293^YFPI152L cells, stably expressing yellow fluorescent protein. Cells were grown for ~ 30 h at 37 °C and ambient CO₂ concentration in biochips mounted into a custom-designed 3D printed carrier and were imaged at a rate of ten images per hour using a fluorescence microscope with environment control infrastructure. Cells in one chip were maintained under static conditions whereas cells in another chip were recurrently perfused with fresh media. Generated image series were quantitatively analyzed using a custom-modified cell detection software. Imaging data averaged from four biological replicates per culturing condition demonstrate that cells cultured under conventional conditions exhibit an exponential growth rate. In contrast, cells cultured in periodic mode exhibited a non-exponential growth rate. Our data clearly indicate differential growth characteristics of cells cultured under periodic versus static conditions highlighting the impact of the culture conditions on the physiology of cells in vitro.

First passage events in biological systems with non-exponential inter-event times

It is often possible to model the dynamics of biological systems as a series of discrete transitions between a finite set of observable states (or compartments). When the residence times in each state, or inter-event times more generally, are exponentially distributed, then one can write a set of ordinary differential equations, which accurately describe the evolution of mean quantities. Non-exponential inter-event times can also be experimentally observed, but are more difficult to analyse mathematically. In this paper, we focus on the computation of first passage events and their probabilities in biological systems with non-exponential inter-event times. We show, with three case studies from Molecular Immunology, Virology and Epidemiology, that significant errors are introduced when drawing conclusions based on the assumption that inter-event times are exponentially distributed. Our approach allows these errors to be avoided with the use of phase-type distributions that approximate arbitrarily distributed inter-event times.

The invasion speed of cell migration models with realistic cell cycle time distributions

Cell proliferation is typically incorporated into stochastic mathematical models of cell migration by assuming that cell divisions occur after an exponentially distributed waiting time. Experimental observations, however, show that this assumption is often far from the real cell cycle time distribution (CCTD). Recent studies have suggested an alternative approach to modelling cell proliferation based on a multi-stage representation of the CCTD. In this paper we investigate the connection between the CCTD and the speed of the collective invasion. We first state a result for a general CCTD, which allows the computation of the invasion speed using the Laplace transform of the CCTD. We use this to deduce the range of speeds for the general case. We then focus on the more realistic case of multi-stage models, using both a stochastic agent-based model and a set of reaction-diffusion equations for the cells' average density. By studying the corresponding travelling wave solutions, we obtain an analytical expression for the speed of invasion for a general N-stage model with identical transition rates, in which case the resulting cell cycle times are Erlang distributed. We show that, for a general N-stage model, the Erlang distribution and the exponential distribution lead to the minimum and maximum invasion speed, respectively. This result allows us to determine the range of possible invasion speeds in terms of the average proliferation time for any multi-stage model.