Patient-derived tumor organoids provide a physiologically relevant 3D disease model for preclinical drug discovery, surpassing the limitations of conventional 2D cell lines. To better capture the dynamic nature of organoid drug responses, we developed a new systematic evaluation method called SCOPE (Systematic Classification of Organoids for Phenotypic Evaluation), harnessing phenotypic assessments from multi-timepoint 3D imaging data. By integrating artificial intelligence (AI)-based image analysis of organoid viability with tracking and mathematical modeling of organoid growth over time, we captured temporal- and dose-dependent dynamics of phenotypic changes, culminating in two novel metrics: a combined growth and viability (GV) score as well as a cytostatic-cytotoxic transition range (CCTR) that separates drug effects on organoid growth and viability. Our approach supports classification of specific drug responses into four distinct phenotypic groups: (1) cytotoxic, (2) cytostatic plus cytotoxic, (3) late cytotoxic, and (4) cytostatic. This novel drug evaluation system can identify previously unknown drug effects or new therapeutic use cases for existing drugs, facilitating the design of alternative therapeutic options to overcome efficacy or drug resistance challenges and improving the clinical applicability of organoid-based drug discovery results.

We consider a supercritical two-type continuous-time linear birth-death process with mutation and selection, in which wild-type individuals give rise to mutant offspring with a larger net growth rate. In this setting, we investigate the “driver” site frequency spectrum (SFS), or the random measure that records mutant allelic frequencies in the population. We examine various regions of the SFS. First, we derive exact moments and prove results concerning the mean behavior of the driver SFS at large times and frequencies. Strong laws of large numbers for the driver SFS are proven by constructing suitable L^2-approximations. These results are extended to the setting in which the fitness increase associated with each clone is random. Next, we allow the frequencies to vary with time to examine the number of “intermediate” and “large” clones. Using this, we find a cutoff frequency at which there are order 1 number of clones. Our results allow for estimation of relevant evolutionary parameters, such as the fitness increase of mutant versus wild-type cells.

Tumor evolution is shaped by cell division, cell death, competition, and constraints imposed by the local microenvironment. Because these dynamics are usually not observed directly, phylogenetic trees inferred from somatic variation in sampled tumor cells can provide an indirect record of the population history that produced the sample. In this paper, we examine whether the distribution of inferred internal branching times exhibits signatures that depend on the underlying tumor growth regime. Specifically, we study the distribution of internal branching times in continuous-time birth-death models of tumor evolution. Exponentially growing populations exhibit a unimodal distribution of internal branching times, with the mode located near the root. In contrast, logistic growth, which models expansion constrained by carrying capacity, yields a substantially more intricate genealogical structure: the distribution of branching times undergoes a systematic transition as the time elapsed since tumor initiation increases. Specifically, this progression shifts from an expansion-dominated phase, through an intermediate early-recent bimodal phase, to a final recent-dominated phase. Extensive simulations of reconstructed tumor genealogies support these theoretical findings.

Therapeutic efficacy of multivalent T cell engagers varies widely across individuals, but the basis for this heterogeneity remains poorly understood. Here, we integrate it in vitro experiments of antitumor immune responses with a mechanistic modeling framework to investigate sources of response variability across T cell donors and TE constructs, focusing on a novel multivalent bispecific T cell engager currently in development. We identify parameter regimes that accurately recapitulate dose-response behaviors across T cell donors and doses, and perform cross-validation studies that demonstrate the model's predictive accuracy. We find that variability in therapy efficacy is governed by the relationship between binding affinity and dose. When dose exceeds the binding affinity, responses are relatively robust across donors; when dose is below the binding affinity, responses are more donor dependent. At smaller doses, the TE-specific shape characteristics of the tumor-binding dose response, its steepness in particular, is a key marker of therapy efficacy. More generally, our integrated modeling and experimental framework offers insights and tools that are applicable to other bispecific T-cell engagers, and provides a quantitative foundation for the systematic, in silico optimization of TE design and dosing strategies. ### Competing Interest Statement The authors have declared no competing interest.

While cancer has traditionally been considered a genetic disease, mounting evidence indicates an important role for non-genetic (epigenetic) mechanisms. Common anti-cancer drugs have recently been observed to induce the adoption of non-genetic drug-tolerant cell states, thereby accelerating the evolution of drug resistance. This confounds conventional high-dose treatment strategies aimed at maximal tumor reduction, since high doses can simultaneously promote non-genetic resistance. In this work, we study optimal dosing of anti-cancer treatment under drug-induced cell plasticity. We show that the optimal dosing strategy steers the tumor to a fixed equilibrium composition between sensitive and tolerant cells, while precisely balancing the trade-off between cell kill and tolerance induction. The optimal equilibrium strategy ranges from applying a low dose continuously to applying the maximum dose intermittently, depending on the dynamics of tolerance induction. We finally discuss how our approach can be integrated with in vitro data to derive patient-specific treatment insights.

Anti-cancer treatment frequently fails due to the evolution of drug resistance. Emerging evidence indicates that in many cases the traditional ‘maximum tolerated dose; treatment paradigm may accelerate resistance evolution. Mathematical modelling can aid in the search for better dosing protocols by enhancing our understanding of the effect of different types of strategies on the evolutionary dynamics of tumour cells and by identifying promising options for pre-clinical and clinical testing. In this chapter, we discuss recent mathematical investigations of the evolution of drug resistance under single-agent continuous and pulsed anti-cancer therapies. We focus on studies involving treatment-induced resistance and phenotypic switching. Our goals are to outline the main ingredients of mathematical models used in recent studies and to extract insights into which model attributes favour continuous versus pulsed schedules and lower versus higher cumulative doses.

Complex interactions between stromal cells, tumor cells and therapies can influence environmental factors that in turn impact anticancer treatment efficacy. Disentangling these phenomena is critical for understanding treatment response and designing effective dosing strategies. We propose a mathematical model for a common tumor-stromal interaction motif where stromal cells secrete factors that promote drug resistance. We demonstrate that the presence of this interaction modulates the therapeutic dose window of efficacy and can lead to nonmonotonic treatment response. We consider combination strategies that target stromal cells and their secretome, and identify strategies that constrain drug concentrations within the efficacious window for long-term response. We explore an experimental dataset from colorectal cancer cells treated with anti-EGFR targeting therapy, cetuximab, where cancer-associated fibroblasts increase epidermal growth factor secretion under treatment. We apply our general approach to identify a critical drug concentration threshold and study effective dosing regimens for single-drug and combination therapies.

PURPOSE:Circulating tumor DNA (ctDNA) assays are promising tools for the prediction of cancer treatment response. Here, we build a framework for the design of ctDNA biomarkers of therapy response that incorporate variations in ctDNA dynamics driven by specific treatment mechanisms. These biomarkers are based on novel proposals for ctDNA sampling protocols, consisting of frequent sampling within a compact time window surrounding therapy initiation-which we hypothesize to hold valuable prognostic information on longer-term treatment response. METHODS:We develop mathematical models of ctDNA kinetics driven by tumor response to several therapy classes and use them to simulate randomized virtual patient cohorts to test candidate biomarkers. RESULTS:Using this approach, we propose specific biomarkers, on the basis of ctDNA longitudinal features, for targeted therapy and radiation therapy. We evaluate and demonstrate the efficacy of these biomarkers in predicting treatment response within a randomized virtual patient cohort data set. CONCLUSION:This study highlights a need for tailoring ctDNA sampling protocols and interpretation methodology to specific biologic mechanisms of therapy response, and it provides a novel modeling and simulation framework for doing so. In addition, it highlights the potential of ctDNA assays for making early, rapid predictions of treatment response within the first days or weeks of treatment and generates hypotheses for further clinical testing.

Multiple myeloma (MM) patients experience repeated cycles of treatment response and relapse, yet despite close monitoring of disease status through M protein measurements, no standard model exists for relapse prediction in MM. We investigate the feasibility of predicting relapse using a hierarchical Bayesian model of subpopulation dynamics by training and testing the model on 229 patients from the IKEMA trial. After observing between 11 and 18 treatment cycles, the model predicted relapse within six cycles with an average sensitivity between 60 and 80 %, and an average specificity between 60 and 90 %. A model of linear extrapolation is preferable when patients have been observed for less than 6 cycles, but for longer observation windows the hierarchical Bayesian model is preferred. Including available baseline and longitudinal covariate information did not improve predictive accuracy. A survival analysis showed that two model parameters separated patients into groups with significantly different PFS ( p < 0.001). Statement of Significance Currently, no standard model exists for relapse prediction in multiple myeloma. A personalized model of M protein development could guide the frequency of follow-up measurements, reduce uncertainty for patients, and give clinicians more time to choose the best subsequent treatment for each patient. Furthermore, models that predict relapse are required to study the effect of changing treatment in advance of relapse rather than in response to it. Our work addresses this need by developing a hierarchical Bayesian model of subpopulation dynamics for prediction of future M protein values. We validate the model on a patient cohort treated with state-of-theart CD38 inhibitor therapy and show that it can accurately predict relapse within the next six treatment cycles, highlighting the promise of mathematical modeling in multiple myeloma and for personalized medicine in general. Declaration of Interests F.S. received honorarium from Sanofi, Janssen, BMS, Oncopeptides, Abbvie, GSK, and Pfizer. The authors declare that they have no other conflicts of interest. ### Competing Interest Statement F.S. received honorarium from Sanofi, Janssen, BMS, Oncopeptides, Abbvie, GSK, and Pfizer. ### Funding Statement E. M. Myklebust, A. Köhn-Luque, and A. Frigessi were supported by the Center for research-based-innovation BigInsight under grant 237718 by the Research Council of Norway. J. Foo and K. Leder were supported by the Fulbright US-Norway Foundation. J. Foo and K. Leder were supported by the University of Oslo-University of Minnesota Norwegian Centennial Chair Grant. J. Foo was supported by the US National Science Foundation under grant number DMS-2052465. K. Leder was supported by the US National Science Foundation under grant number CMMI-2228034. We acknowledge funding from the Research Council of Norway through projects DL: Pipeline for individually tailoring new treatments in hematological cancers (PINpOINT) under project number 294916, and INTPART-International Partnerships for Excellent Education and Research under project number 309273. The authors also acknowledge the Centre for Digital Life Norway for supporting the partner project PINpOINT. ### Author Declarations I confirm all relevant ethical guidelines have been followed, and any necessary IRB and/or ethics committee approvals have been obtained. Yes The details of the IRB/oversight body that provided approval or exemption for the research described are given below: The study protocol of the IKEMA trial ([NCT03275285][1]) was approved by the Institutional Ethics Committee or independent review board for each center. I confirm that all necessary patient/participant consent has been obtained and the appropriate institutional forms have been archived, and that any patient/participant/sample identifiers included were not known to anyone (e.g., hospital staff, patients or participants themselves) outside the research group so cannot be used to identify individuals. Yes I understand that all clinical trials and any other prospective interventional studies must be registered with an ICMJE-approved registry, such as ClinicalTrials.gov. I confirm that any such study reported in the manuscript has been registered and the trial registration ID is provided (note: if posting a prospective study registered retrospectively, please provide a statement in the trial ID field explaining why the study was not registered in advance). Yes I have followed all appropriate research reporting guidelines, such as any relevant EQUATOR Network research reporting checklist(s) and other pertinent material, if applicable. Yes The code for this project is available at https://github.com/evenmm/mm-predict-ikema. Data from the IKEMA trial ([NCT03275285][1]) can be requested through the data-sharing platform Vivli. <https://search.vivli.org/studyDetails/fromSearch/42d32e57-51ce-4d54-9613-eb3f7c73d30e> <https://github.com/evenmm/mm-predict-ikema> [1]: /lookup/external-ref?link_type=CLINTRIALGOV&access_num=NCT03275285&atom=%2Fmedrxiv%2Fearly%2F2024%2F05%2F06%2F2024.05.02.24306607.atom

Patient-derived tumor organoids (PDTOs) are novel cellular models that maintain the genetic, phenotypic and structural features of patient tumor tissue and are useful for studying tumorigenesis and drug response. When integrated with advanced 3D imaging and analysis techniques, PDTOs can be used to establish physiologically relevant high-throughput and high-content drug screening platforms that support the development of patient-specific treatment strategies. However, in order to effectively leverage high-throughput PDTO observations for clinical predictions, it is critical to establish a quantitative understanding of the basic properties and variability of organoid growth dynamics. In this work, we introduced an innovative workflow for analyzing and understanding PDTO growth dynamics, by integrating a high-throughput imaging deep learning platform with mathematical modeling, incorporating flexible growth laws and variable dormancy times. We applied the workflow to colon cancer organoids and demonstrated that organoid growth is well-described by the Gompertz model of growth. Our analysis showed significant intrapatient heterogeneity in PDTO growth dynamics, with the initial exponential growth rate of an organoid following a lognormal distribution within each dataset. The level of intrapatient heterogeneity varied between patients, as did organoid growth rates and dormancy times of single seeded cells. Our work contributes to an emerging understanding of the basic growth characteristics of PDTOs, and it highlights the heterogeneity in organoid growth both within and between patients. These results pave the way for further modeling efforts aimed at predicting treatment response dynamics and drug resistance timing.

Circulating tumor DNA assays are promising tools for the prediction of cancer treatment response. Here, we build a framework for the design of ctDNA biomarkers of therapy response that incorporate variations in ctDNA dynamics driven by specific treatment mechanisms. We develop mathematical models of ctDNA kinetics driven by tumor response to several therapy classes, and utilize them to simulate randomized virtual patient cohorts to test candidate biomarkers. Using this approach, we propose specific biomarkers, based on ctDNA longitudinal features, for targeted therapy, chemotherapy and radiation therapy. We evaluate and demonstrate the efficacy of these biomarkers in predicting treatment response within a randomized virtual patient cohort dataset. These biomarkers are based on novel proposals for ctDNA sampling protocols, consisting of frequent sampling within a compact time window surrounding therapy initiation - which we hypothesize to hold valuable prognostic information on longer-term treatment response. This study highlights a need for tailoring ctDNA sampling protocols and interpretation methodology to specific biological mechanisms of therapy response, and it provides a novel modeling and simulation framework for doing so. In addition, it highlights the potential of ctDNA assays for making early, rapid predictions of treatment response within the first days or weeks of treatment, and generates hypotheses for further clinical testing.

Tumor heterogeneity is a complex and widely recognized trait that poses significant challenges in developing effective cancer therapies. In particular, many tumors harbor a variety of subpopulations with distinct therapeutic response characteristics. Characterizing this heterogeneity by determining the subpopulation structure within a tumor enables more precise and successful treatment strategies. In our prior work, we developed PhenoPop, a computational framework for unravelling the drug-response subpopulation structure within a tumor from bulk high-throughput drug screening data. However, the deterministic nature of the underlying models driving PhenoPop restricts the model fit and the information it can extract from the data. As an advancement, we propose a stochastic model based on the linear birth-death process to address this limitation. Our model can formulate a dynamic variance along the horizon of the experiment so that the model uses more information from the data to provide a more robust estimation. In addition, the newly proposed model can be readily adapted to situations where the experimental data exhibits a positive time correlation. We test our model on simulated data (in silico) and experimental data (in vitro), which supports our argument about its advantages.

Predicting cancer dynamics under treatment is challenging due to high inter-patient heterogeneity, lack of predictive biomarkers, and sparse and noisy longitudinal data. Mathematical models can summarize cancer dynamics by a few interpretable parameters per patient. Machine learning methods can then be trained to predict the model parameters from baseline covariates, but do not account for uncertainty in the parameter estimates. Instead, hierarchical Bayesian modeling can model the relationship between baseline covariates to longitudinal measurements via mechanistic parameters while accounting for uncertainty in every part of the model. The mapping from baseline covariates to model parameters can be modeled in several ways. A linear mapping simplifies inference but fails to capture nonlinear covariate effects and scale poorly for interaction modeling when the number of covariates is large. In contrast, Bayesian neural networks can potentially discover interactions between covariates automatically, but at a substantial cost in computational complexity. In this work, we develop a hierarchical Bayesian model of subpopulation dynamics that uses baseline covariate information to predict cancer dynamics under treatment, inspired by cancer dynamics in multiple myeloma (MM), where serum M protein is a well-known proxy of tumor burden. As a working example, we apply the model to a simulated dataset and compare its ability to predict M protein trajectories to a model with linear covariate effects. Our results show that the Bayesian neural network covariate effect model predicts cancer dynamics more accurately than a linear covariate effect model when covariate interactions are present. The framework can also be applied to other types of cancer or other time series prediction problems that can be described with a parametric model.

Tumor recurrence, driven by the evolution of drug resistance is a major barrier to therapeutic success in cancer. Resistance is often caused by genetic alterations such as point mutation, which refers to the modification of a single genomic base pair, or gene amplification, which refers to the duplication of a region of DNA that contains a gene. Here we investigate the dependence of tumor recurrence dynamics on these mechanisms of resistance, using stochastic multi-type branching process models. We derive tumor extinction probabilities and deterministic estimates for the tumor recurrence time, defined as the time when an initially drug sensitive tumor surpasses its original size after developing resistance. For models of amplification-driven and mutation-driven resistance, we prove law of large numbers results regarding the convergence of the stochastic recurrence times to their mean. Additionally, we prove sufficient and necessary conditions for a tumor to escape extinction under the gene amplification model, discuss behavior under biologically relevant parameters, and compare the recurrence time and tumor composition in the mutation and amplification models both analytically and using simulations. In comparing these mechanisms, we find that the ratio between recurrence times driven by amplification vs. mutation depends linearly on the number of amplification events required to acquire the same degree of resistance as a mutation event, and we find that the relative frequency of amplification and mutation events plays a key role in determining the mechanism under which recurrence is more rapid. In the amplification-driven resistance model, we also observe that increasing drug concentration leads to a stronger initial reduction in tumor burden, but that the eventual recurrent tumor population is less heterogeneous, more aggressive, and harbors higher levels of drug-resistance.

Over 80% of human cancers originate from the epithelium, which covers the outer and inner surfaces of organs and blood vessels. In stratified epithe-lium, the bottom layers are occupied by stem and stem-like cells that contin-ually divide and replenish the upper layers. In this work, we study the spread of premalignant mutant clones and cancer initiation in stratified epithelium, using the biased voter model on stacked two-dimensional lattices. Our main result is an estimate of the propagation speed of a premalignant mutant clone, which is asymptotically precise in the cancer-relevant weak-selection limit. We use our main result to study cancer initiation under a two-step mutational model of cancer, which includes computing the distributions of the time of cancer initiation and the size of the premalignant clone giving rise to cancer. Our work quantifies the effect of epithelial tissue thickness on the process of carcinogenesis, thereby contributing to an emerging understanding of the spatial evolutionary dynamics of cancer.

Recent evidence suggests that nongenetic (epigenetic) mechanisms play an important role at all stages of cancer evolution. In many cancers, these mechanisms have been observed to induce dynamic switching between two or more cell states, which commonly show differential responses to drug treatments. To understand how these cancers evolve over time, and how they respond to treatment, we need to understand the state-dependent rates of cell proliferation and phenotypic switching. In this work, we propose a rigorous statistical framework for estimating these parameters, using data from commonly performed cell line experiments, where phenotypes are sorted and expanded in culture. The framework explicitly models the stochastic dynamics of cell division, cell death and phenotypic switching, and it provides likelihood-based confidence intervals for the model parameters. The input data can be either the fraction of cells or the number of cells in each state at one or more time points. Through a combination of theoretical analysis and numerical simulations, we show that when cell fraction data is used, the rates of switching may be the only parameters that can be estimated accurately. On the other hand, using cell number data enables accurate estimation of the net division rate for each phenotype, and it can even enable estimation of the state-dependent rates of cell division and cell death. We conclude by applying our framework to a publicly available dataset.

Tumor heterogeneity is an important driver of treatment failure in cancer since therapies often select for drug-tolerant or drug-resistant cellular subpopulations that drive tumor growth and recurrence. Profiling the drug-response heterogeneity of tumor samples using traditional genomic deconvolution methods has yielded limited results, due in part to the imperfect mapping between genomic variation and functional characteristics. Here, we leverage mechanistic population modeling to develop a statistical framework for profiling phenotypic heterogeneity from standard drug-screen data on bulk tumor samples. This method, called PhenoPop, reliably identifies tumor subpopulations exhibiting differential drug responses and estimates their drug sensitivities and frequencies within the bulk population. We apply PhenoPop to synthetically generated cell populations, mixed cell-line experiments, and multiple myeloma patient samples and demonstrate how it can provide individualized predictions of tumor growth under candidate therapies. This methodology can also be applied to deconvolution problems in a variety of biological settings beyond cancer drug response.

This file contains four sections with mathematical and methodological details, including one supplemental figure. Section S1: Based on the mesoscopic mathematical model, the likelihood function for Bayesian parameter inference is derived. Section S2: A priori parameter estimates are justified based on experimental findings from the literature. Section S3: The probabilistic distributions of field size and multiplicity are derived in detail. Section S4: Age-specific incidence data is compared to the classical Armitage-Doll model of multistage carcinogenesis. Figure S1: The Armitage-Doll model fit to age-specific incidence data is visualized

The spread of an advantageous mutation through a population is of fundamental interest in population genetics. While the classical Moran model is formulated for a well-mixed population, it has long been recognized that in real-world applications, the population usually has an explicit spatial structure which can significantly influence the dynamics. In the context of cancer initiation in epithelial tissue, several recent works have analyzed the dynamics of advantageous mutant spread on integer lattices, using the biased voter model from particle systems theory. In this spatial version of the Moran model, individuals first reproduce according to their fitness and then replace a neighboring individual. From a biological standpoint, the opposite dynamics, where individuals first die and are then replaced by a neighboring individual according to its fitness, are equally relevant. Here, we investigate this death-birth analogue of the biased voter model. We construct the process mathematically, derive the associated dual process, establish bounds on the survival probability of a single mutant, and prove that the process has an asymptotic shape. We also briefly discuss alternative birth-death and death-birth dynamics, depending on how the mutant fitness advantage affects the dynamics. We show that birth-death and death-birth formulations of the biased voter model are equivalent when fitness affects the former event of each update of the model, whereas the birth-death model is fundamentally different from the death-birth model when fitness affects the latter event.