Hybrid Automata Library: A flexible platform for hybrid modeling with real-time visualization.

The Hybrid Automata Library (HAL) is a Java Library developed for use in mathematical oncology modeling. It is made of simple, efficient, generic components that can be used to model complex spatial systems. HAL's components can broadly be classified into: on- and off-lattice agent containers, finite difference diffusion fields, a GUI building system, and additional tools and utilities for computation and data collection. These components are designed to operate independently and are standardized to make them easy to interface with one another. As a demonstration of how modeling can be simplified using our approach, we have included a complete example of a hybrid model (a spatial model with interacting agent-based and PDE components). HAL is a useful asset for researchers who wish to build performant 1D, 2D and 3D hybrid models in Java, while not starting entirely from scratch. It is available on GitHub at https://github.com/MathOnco/HAL under the MIT License. HAL requires the Java JDK version 1.8 or later to compile and run the source code.

From cells to tissue: How cell scale heterogeneity impacts glioblastoma growth and treatment response.

Glioblastomas are aggressive primary brain tumors known for their inter- and intratumor heterogeneity. This disease is uniformly fatal, with intratumor heterogeneity the major reason for treatment failure and recurrence. Just like the nature vs nurture debate, heterogeneity can arise from intrinsic or environmental influences. Whilst it is impossible to clinically separate observed behavior of cells from their environmental context, using a mathematical framework combined with multiscale data gives us insight into the relative roles of variation from different sources. To better understand the implications of intratumor heterogeneity on therapeutic outcomes, we created a hybrid agent-based mathematical model that captures both the overall tumor kinetics and the individual cellular behavior. We track single cells as agents, cell density on a coarser scale, and growth factor diffusion and dynamics on a finer scale over time and space. Our model parameters were fit utilizing serial MRI imaging and cell tracking data from ex vivo tissue slices acquired from a growth-factor driven glioblastoma murine model. When fitting our model to serial imaging only, there was a spectrum of equally-good parameter fits corresponding to a wide range of phenotypic behaviors. When fitting our model using imaging and cell scale data, we determined that environmental heterogeneity alone is insufficient to match the single cell data, and intrinsic heterogeneity is required to fully capture the migration behavior. The wide spectrum of in silico tumors also had a wide variety of responses to an application of an anti-proliferative treatment. Recurrent tumors were generally less proliferative than pre-treatment tumors as measured via the model simulations and validated from human GBM patient histology. Further, we found that all tumors continued to grow with an anti-migratory treatment alone, but the anti-proliferative/anti-migratory combination generally showed improvement over an anti-proliferative treatment alone. Together our results emphasize the need to better understand the underlying phenotypes and tumor heterogeneity present in a tumor when designing therapeutic regimens.

Lesion Dynamics Under Varying Paracrine PDGF Signaling in Brain Tissue.

Paracrine PDGF signaling is involved in many processes in the body, both normal and pathological, including embryonic development, angiogenesis, and wound healing as well as liver fibrosis, atherosclerosis, and cancers. We explored this seemingly dual (normal and pathological) role of PDGF mathematically by modeling the release of PDGF in brain tissue and then varying the dynamics of this release. Resulting simulations show that by varying the dynamics of a PDGF source, our model predicts three possible outcomes for PDGF-driven cellular recruitment and lesion growth: (1) localized, short duration of growth, (2) localized, chronic growth, and (3) widespread chronic growth. Further, our model predicts that the type of response is much more sensitive to the duration of PDGF exposure than the maximum level of that exposure. This suggests that extended duration of paracrine PDGF signal during otherwise normal processes could potentially lead to lesions having a phenotype consistent with pathologic conditions.

Methods for determining key components in a mathematical model for tumor-immune dynamics in multiple myeloma.

In this work, we analyze a mathematical model we introduced previously for the dynamics of multiple myeloma and the immune system. We focus on four main aspects: (1) obtaining and justifying ranges and values for all parameters in the model; (2) determining a subset of parameters to which the model is most sensitive; (3) determining which parameters in this subset can be uniquely estimated given certain types of data; and (4) exploring the model numerically. Using global sensitivity analysis techniques, we found that the model is most sensitive to certain growth, loss, and efficacy parameters. This analysis provides the foundation for a future application of the model: prediction of optimal combination regimens in patients with multiple myeloma.

Simulating PDGF-Driven Glioma Growth and Invasion in an Anatomically Accurate Brain Domain.

Gliomas are the most common of all primary brain tumors. They are characterized by their diffuse infiltration of the brain tissue and are uniformly fatal, with glioblastoma being the most aggressive form of the disease. In recent years, the over-expression of platelet-derived growth factor (PDGF) has been shown to produce tumors in experimental rodent models that closely resemble this human disease, specifically the proneural subtype of glioblastoma. We have previously modeled this system, focusing on the key attribute of these experimental tumors-the "recruitment" of oligodendroglial progenitor cells (OPCs) to participate in tumor formation by PDGF-expressing retrovirally transduced cells-in one dimension, with spherical symmetry. However, it has been observed that these recruitable progenitor cells are not uniformly distributed throughout the brain and that tumor cells migrate at different rates depending on the material properties in different regions of the brain. Here we model the differential diffusion of PDGF-expressing and recruited cell populations via a system of partial differential equations with spatially variable diffusion coefficients and solve the equations in two spatial dimensions on a mouse brain atlas using a flux-differencing numerical approach. Simulations of our in silico model demonstrate qualitative agreement with the observed tumor distribution in the experimental animal system. Additionally, we show that while there are higher concentrations of OPCs in white matter, the level of recruitment of these plays little role in the appearance of "white matter disease," where the tumor shows a preponderance for white matter. Instead, simulations show that this is largely driven by the ratio of the diffusion rate in white matter as compared to gray. However, this ratio has less effect on the speed of tumor growth than does the degree of OPC recruitment in the tumor. It was observed that tumor simulations with greater degrees of recruitment grow faster and develop more nodular tumors than if there is no recruitment at all, similar to our prior results from implementing our model in one dimension. Combined, these results show that recruitment remains an important consideration in understanding and slowing glioma growth.

To modulate or to skip: De-escalating PARP inhibitor maintenance therapy in ovarian cancer using adaptive therapy.

Toxicity and emerging drug resistance pose important challenges in poly-adenosine ribose polymerase inhibitor (PARPi) maintenance therapy of ovarian cancer. We propose that adaptive therapy, which dynamically reduces treatment based on the tumor dynamics, might alleviate both issues. Utilizing in vitro time-lapse microscopy and stepwise model selection, we calibrate and validate a differential equation mathematical model, which we leverage to test different plausible adaptive treatment schedules. Our model indicates that adjusting the dosage, rather than skipping treatments, is more effective at reducing drug use while maintaining efficacy due to a delay in cell kill and a diminishing dose-response relationship. In vivo pilot experiments confirm this conclusion. Although our focus is toxicity mitigation, reducing drug use may also delay resistance. This study enhances our understanding of PARPi treatment scheduling and illustrates the first steps in developing adaptive therapies for new treatment settings. A record of this paper's transparent peer review process is included in the supplemental information.

Adaptive therapy for ovarian cancer: An integrated approach to PARP inhibitor scheduling.

Toxicity and emerging drug resistance are important challenges in PARP inhibitor (PARPi) treatment of ovarian cancer. Recent research has shown that evolutionary-inspired treatment algorithms which adapt treatment to the tumor's treatment response (adaptive therapy) can help to mitigate both. Here, we present a first step in developing an adaptive therapy protocol for PARPi treatment by combining mathematical modelling and wet-lab experiments to characterize the cell population dynamics under different PARPi schedules. Using data from in vitro Incucyte Zoom time-lapse microscopy experiments and a step-wise model selection process we derive a calibrated and validated ordinary differential equation model, which we then use to test different plausible adaptive treatment schedules. Our model can accurately predict the in vitro treatment dynamics, even to new schedules, and suggests that treatment modifications need to be carefully timed, or one risks losing control over tumour growth, even in the absence of any resistance. This is because our model predicts that multiple rounds of cell division are required for cells to acquire sufficient DNA damage to induce apoptosis. As a result, adaptive therapy algorithms that modulate treatment but never completely withdraw it are predicted to perform better in this setting than strategies based on treatment interruptions. Pilot experiments in vivo confirm this conclusion. Overall, this study contributes to a better understanding of the impact of scheduling on treatment outcome for PARPis and showcases some of the challenges involved in developing adaptive therapies for new treatment settings.

Intermetastatic and Intrametastatic Heterogeneity Shapes Adaptive Therapy Cycling Dynamics.

Adaptive therapies that alternate between drug applications and drug-free vacations can exploit competition between sensitive and resistant cells to maximize the time to progression. However, optimal dosing schedules depend on the properties of metastases, which are often not directly measurable in clinical practice. Here, we proposed a framework for estimating features of metastases through tumor response dynamics during the first adaptive therapy treatment cycle. Longitudinal prostate-specific antigen (PSA) levels in 16 patients with metastatic castration-resistant prostate cancer undergoing adaptive androgen deprivation treatment were analyzed to investigate relationships between cycle dynamics and clinical variables such as Gleason score, the change in the number of metastases over a cycle, and the total number of cycles over the course of treatment. The first cycle of adaptive therapy, which consists of a response period (applying therapy until 50% PSA reduction), and a regrowth period (removing treatment until reaching initial PSA levels), delineated several features of the computational metastatic system: larger metastases had longer cycles; a higher proportion of drug-resistant cells slowed the cycles; and a faster cell turnover rate sped up drug response time and slowed regrowth time. The number of metastases did not affect cycle times, as response dynamics were dominated by the largest tumors rather than the aggregate. In addition, systems with higher intermetastasis heterogeneity responded better to continuous therapy and correlated with dynamics from patients with high or low Gleason scores. Conversely, systems with higher intrametastasis heterogeneity responded better to adaptive therapy and correlated with dynamics from patients with intermediate Gleason scores. SIGNIFICANCE: Multiscale mathematical modeling combined with biomarker dynamics during adaptive therapy helps identify underlying features of metastatic cancer to inform treatment decisions.

A survey of open questions in adaptive therapy: Bridging mathematics and clinical translation.

Adaptive therapy is a dynamic cancer treatment protocol that updates (or 'adapts') treatment decisions in anticipation of evolving tumor dynamics. This broad term encompasses many possible dynamic treatment protocols of patient-specific dose modulation or dose timing. Adaptive therapy maintains high levels of tumor burden to benefit from the competitive suppression of treatment-sensitive subpopulations on treatment-resistant subpopulations. This evolution-based approach to cancer treatment has been integrated into several ongoing or planned clinical trials, including treatment of metastatic castrate resistant prostate cancer, ovarian cancer, and BRAF-mutant melanoma. In the previous few decades, experimental and clinical investigation of adaptive therapy has progressed synergistically with mathematical and computational modeling. In this work, we discuss 11 open questions in cancer adaptive therapy mathematical modeling. The questions are split into three sections: (1) integrating the appropriate components into mathematical models (2) design and validation of dosing protocols, and (3) challenges and opportunities in clinical translation.

Spatial structure impacts adaptive therapy by shaping intra-tumoral competition.

Background: Adaptive therapy aims to tackle cancer drug resistance by leveraging resource competition between drug-sensitive and resistant cells. Here, we present a theoretical study of intra-tumoral competition during adaptive therapy, to investigate under which circumstances it will be superior to aggressive treatment. Methods: We develop and analyse a simple, 2-D, on-lattice, agent-based tumour model in which cells are classified as fully drug-sensitive or resistant. Subsequently, we compare this model to its corresponding non-spatial ordinary differential equation model, and fit it to longitudinal prostate-specific antigen data from 65 prostate cancer patients undergoing intermittent androgen deprivation therapy following biochemical recurrence. Results: Leveraging the individual-based nature of our model, we explicitly demonstrate competitive suppression of resistance during adaptive therapy, and examine how different factors, such as the initial resistance fraction or resistance costs, alter competition. This not only corroborates our theoretical understanding of adaptive therapy, but also reveals that competition of resistant cells with each other may play a more important role in adaptive therapy in solid tumours than was previously thought. To conclude, we present two case studies, which demonstrate the implications of our work for: (i) mathematical modelling of adaptive therapy, and (ii) the intra-tumoral dynamics in prostate cancer patients during intermittent androgen deprivation treatment, a precursor of adaptive therapy. Conclusion: Our work shows that the tumour's spatial architecture is an important factor in adaptive therapy and provides insights into how adaptive therapy leverages both inter- and intra-specific competition to control resistance.

A Phase 1b Adaptive Androgen Deprivation Therapy Trial in Metastatic Castration Sensitive Prostate Cancer.

Background: We hypothesize that cancer survival can be improved through adapting treatment strategies to cancer evolutionary dynamics and conducted a phase 1b study in metastatic castration sensitive prostate cancer (mCSPC). Methods: Men with asymptomatic mCSPC were enrolled and proceeded with a treatment break after achieving > 75% PSA decline with LHRH analog plus an NHA. ADT was restarted at the time of PSA or radiographic progression and held again after achieving >50% PSA decline. This on-off cycling of ADT continued until on treatment imaging progression. Results: At data cut off in August 2022, only 2 of the 16 evaluable patients were off study due to imaging progression at 28 months from first dose of LHRH analog for mCSPC. Two additional patients showed PSA progression at 12.4 and 20.5 months and remain on trial. Since none of the 16 patients developed imaging progression at 12 months, the study succeeded in its primary objective of feasibility. The secondary endpoints of median time to PSA progression and median time to radiographic progression have not been reached at a median follow up of 26 months. Conclusions: It is feasible to use an individual's PSA response and testosterone levels to guide intermittent ADT in mCSPC.

Cycling hypoxia selects for constitutive HIF stabilization.

Tumors experience temporal and spatial fluctuations in oxygenation. Hypoxia inducible transcription factors (HIF-α) respond to low levels of oxygen and induce re-supply oxygen. HIF-α stabilization is typically facultative, induced by hypoxia and reduced by normoxia. In some cancers, HIF-α stabilization becomes constitutive under normoxia. We develop a mathematical model that predicts how fluctuating oxygenation affects HIF-α stabilization and impacts net cell proliferation by balancing the base growth rate, the proliferative cost of HIF-α expression, and the mortality from not expressing HIF-α during hypoxia. We compare optimal net cell proliferation rate between facultative and constitutive HIF-α regulation in environments with different oxygen profiles. We find that that facultative HIF-α regulation promotes greater net cell proliferation than constitutive regulation with stochastic or slow periodicity in oxygenation. However, cell fitness is nearly identical for both HIF-α regulation strategies under rapid periodic oxygenation fluctuations. The model thus indicates that cells constitutively expressing HIF-α may be at a selective advantage when the cost of expression is low. In cancer, this condition is known as pseudohypoxia or the "Warburg Effect". We conclude that rapid and regular cycling of oxygenation levels selects for pseudohypoxia, and that this is consistent with the ecological theory of optimal defense.

The impact of proliferation-migration tradeoffs on phenotypic evolution in cancer.

Tumors are not static masses of cells but dynamic ecosystems where cancer cells experience constant turnover and evolve fitness-enhancing phenotypes. Selection for different phenotypes may vary with (1) the tumor niche (edge or core), (2) cell turnover rates, (3) the nature of the tradeoff between traits, and (4) whether deaths occur in response to demographic or environmental stochasticity. Using a spatially-explicit agent-based model, we observe how two traits (proliferation rate and migration speed) evolve under different tradeoff conditions with different turnover rates. Migration rate is favored over proliferation at the tumor's edge and vice-versa for the interior. Increasing cell turnover rates slightly slows tumor growth but accelerates the rate of evolution for both proliferation and migration. The absence of a tradeoff favors ever higher values for proliferation and migration, while a convex tradeoff tends to favor proliferation, often promoting the coexistence of a generalist and specialist phenotype. A concave tradeoff favors migration at low death rates, but switches to proliferation at higher death rates. Mortality via demographic stochasticity favors proliferation, and environmental stochasticity favors migration. While all of these diverse factors contribute to the ecology, heterogeneity, and evolution of a tumor, their effects may be predictable and empirically accessible.

Spatial Heterogeneity and Evolutionary Dynamics Modulate Time to Recurrence in Continuous and Adaptive Cancer Therapies.

Treatment of advanced cancers has benefited from new agents that supplement or bypass conventional therapies. However, even effective therapies fail as cancer cells deploy a wide range of resistance strategies. We propose that evolutionary dynamics ultimately determine survival and proliferation of resistant cells. Therefore, evolutionary strategies should be used with conventional therapies to delay or prevent resistance. Using an agent-based framework to model spatial competition among sensitive and resistant populations, we applied antiproliferative drug treatments to varying ratios of sensitive and resistant cells. We compared a continuous maximum-tolerated dose schedule with an adaptive schedule aimed at tumor control via competition between sensitive and resistant cells. Continuous treatment cured mostly sensitive tumors, but with any resistant cells, recurrence was inevitable. We identified two adaptive strategies that control heterogeneous tumors: dose modulation controls most tumors with less drug, while a more vacation-oriented schedule can control more invasive tumors. These findings offer potential modifications to treatment regimens that may improve outcomes and reduce resistance and recurrence.Significance: By using drug dose modulation or treatment vacations, adaptive therapy strategies control the emergence of tumor drug resistance by spatially suppressing less fit resistant populations in favor of treatment sensitive ones. Cancer Res; 78(8); 2127-39. ©2018 AACR.

Comparative dynamics of microglial and glioma cell motility at the infiltrative margin of brain tumours.

Microglia are a major cellular component of gliomas, and abundant in the centre of the tumour and at the infiltrative margins. While glioma is a notoriously infiltrative disease, the dynamics of microglia and glioma migratory patterns have not been well characterized. To investigate the migratory behaviour of microglia and glioma cells at the infiltrative edge, we performed two-colour time-lapse fluorescence microscopy of brain slices generated from a platelet-derived growth factor-B (PDGFB)-driven rat model of glioma, in which glioma cells and microglia were each labelled with one of two different fluorescent markers. We used mathematical techniques to analyse glioma cells and microglia motility with both single cell tracking and particle image velocimetry (PIV). Our results show microglia motility is strongly correlated with the presence of glioma, while the correlation of the speeds of glioma cells and microglia was variable and weak. Additionally, we showed that microglia and glioma cells exhibit different types of diffusive migratory behaviour. Microglia movement fit a simple random walk, while glioma cell movement fits a super diffusion pattern. These results show that glioma cells stimulate microglia motility at the infiltrative margins, creating a correlation between the spatial distribution of glioma cells and the pattern of microglia motility.

Bridging population and tissue scale tumor dynamics: a new paradigm for understanding differences in tumor growth and metastatic disease.

To provide a better understanding of the relationship between primary tumor growth rates and metastatic burden, we present a method that bridges tumor growth dynamics at the population level, extracted from the SEER database, to those at the tissue level. Specifically, with this method, we are able to relate estimates of tumor growth rates and metastatic burden derived from a population-level model to estimates of the primary tumor vascular response and the circulating tumor cell (CTC) fraction derived from a tissue-level model. Variation in the population-level model parameters produces differences in cancer-specific survival and cure fraction. Variation in the tissue-level model parameters produces different primary tumor dynamics that subsequently lead to different growth dynamics of the CTCs. Our method to bridge the population and tissue scales was applied to lung and breast cancer separately, and the results were compared. The population model suggests that lung tumors grow faster and shed a significant number of lethal metastatic cells at small sizes, whereas breast tumors grow slower and do not significantly shed lethal metastatic cells until becoming larger. Although the tissue-level model does not explicitly model the metastatic population, we are able to disengage the direct dependency of the metastatic burden on primary tumor growth by introducing the CTC population as an intermediary and assuming dependency. We calibrate the tissue-level model to produce results consistent with the population model while also revealing a more dynamic relationship between the primary tumor and the CTCs. This leads to exponential tumor growth in lung and power law tumor growth in breast. We conclude that the vascular response of the primary tumor is a major player in the dynamics of both the primary tumor and the CTCs, and is significantly different in breast and lung cancer.

Improving treatment strategies for patients with metastatic castrate resistant prostate cancer through personalized computational modeling.

Metastatic castrate resistant prostate cancer (mCRPC) is responsible for the majority of prostate cancer deaths with the median survival after diagnosis being 2 years. The metastatic lesions often arise in the skeleton, and current treatment options are primarily palliative. Using guidelines set forth by the National Comprehensive Cancer Network (NCCN), the medical oncologist has a number of choices available to treat the metastases. However, the sequence of those treatments is largely dependent on the patient history, treatment response and preferences. We posit that the utilization of personalized computational models and treatment optimization algorithms based on patient specific parameters could significantly enhance the oncologist's ability to choose an optimized sequence of available therapies to maximize overall survival. In this perspective, we used an integrated team approach involving clinicians, researchers, and mathematicians, to generate an example of how computational models and genetic algorithms can be utilized to predict the response of heterogeneous mCRPCs in bone to varying sequences of standard and targeted therapies. The refinement and evolution of these powerful models will be critical for extending the overall survival of men diagnosed with mCRPC.

Evolution of intratumoral phenotypic heterogeneity: the role of trait inheritance.

A tumour is a heterogeneous population of cells that competes for limited resources. In the clinic, we typically probe the tumour by biopsy, and then characterize it by the dominant genetic clone. But genotypes are only the first link in the chain of hierarchical events that leads to a specific cell phenotype. The relationship between genotype and phenotype is not simple, and the so-called genotype to phenotype map is poorly understood. Many genotypes can produce the same phenotype, so genetic heterogeneity may not translate directly to phenotypic heterogeneity. We therefore choose to focus on the functional endpoint, the phenotype as defined by a collection of cellular traits (e.g. proliferative and migratory ability). Here, we will examine how phenotypic heterogeneity evolves in space and time and how the way in which phenotypes are inherited will drive this evolution. A tumour can be thought of as an ecosystem, which critically means that we cannot just consider it as a collection of mutated cells but more as a complex system of many interacting cellular and microenvironmental elements. At its simplest, a growing tumour with increased proliferation capacity must compete for space as a limited resource. Hypercellularity leads to a contact-inhibited core with a competitive proliferating rim. Evolution and selection occurs, and an individual cell's capacity to survive and propagate is determined by its combination of traits and interaction with the environment. With heterogeneity in phenotypes, the clone that will dominate is not always obvious as there are both local interactions and global pressures. Several combinations of phenotypes can coexist, changing the fitness of the whole. To understand some aspects of heterogeneity in a growing tumour, we build an off-lattice agent-based model consisting of individual cells with assigned trait values for proliferation and migration rates. We represent heterogeneity in these traits with frequency distributions and combinations of traits with density maps. How the distributions change over time is dependent on how traits are passed on to progeny cells, which is our main enquiry. We bypass the translation of genetics to behaviour by focusing on the functional end result of inheritance of the phenotype combined with the environmental influence of limited space.

First order phase transition and hysteresis in a cell's maintenance of the membrane potential--An essential role for the inward potassium rectifiers.

Hysteretic behavior is found experimentally in the transmembrane potential at low extracellular potassium in mouse lumbrical muscle cells. Adding isoprenaline to the external medium eliminates the bistable, hysteretic region. The system can be modeled mathematically and understood analytically with and without isoprenaline. Inward rectifying potassium channels appear to be essential for the bistability. Relations are derived to express the dimensions of the bistable area in terms of system parameters. The selective advantage and evolutionary origin of inward rectifying channels and hysteretic behavior is discussed.

Ion-channel-like behavior in lipid bilayer membranes at the melting transition.

It is well known that at the gel-liquid phase-transition temperature a lipid bilayer membrane exhibits an increased ion permeability. We analyze the quantized currents in which the increased permeability presents itself. The open time histogram shows a "-3/2" power law which implies an open-closed transition rate that decreases like k(t)∝t(-1) as time evolves. We propose a "pore freezing" model to explain the observations. We discuss how this model also leads to the 1/fα noise that is commonly observed in currents across biological and artificial membranes.