Blood vessel adaptation plays an important role in maintaining healthy and well-oxygenated tissue throughout the body. This is especially true for the brain. To better characterize how blood flow changes when the brain suffers a major arterial occlusion (e.g. during a stroke) and to identify major factors that may affect flow restoration to downstream regions, we created a mathematical model of blood flow in the brain. The network is modeled as multiple larger vessels interconnected with multiple compartments of smaller vessels with each compartment consisting of identical vessels situated in parallel. The model further includes vessel adaptation in response to changes in pressure (myogenic response), wall shear stress (shear response), and oxygen saturation (metabolic response). By varying tissue oxygen demand and incoming pressure, we are able to identify that the number of collateral vessels moving flow from unobstructed to obstructed regions is a major factor. We also predicted a loss of normal function particularly reflected by a shift in the “autoregulation curve”, a curve that reflects the ability of vessels to reasonably respond to increases in pressure. Such results were consistent with experiment and reinforce the appropriateness of treatments that raise flow and oxygenation by minimizing tissue oxygen demand and raising vascular pressure.

Blood filtration occurs in renal capillary tufts called glomeruli, the internal structure of which bears on questions of function, development, evolution, and pathology. Due to the low resolution and labor-intensity of imaging technology, only a handful of studies reaching back decades have examined the spatial structure of glomerular capillaries. Several common features have been described, including lobular topology, plausibly associated with robustness to vascular damage, and circuitous geometry, hypothesized to ensure consistent filtration. However, these properties have been neither mathematically defined nor statistically confirmed. Recent developments in serial scanning electron microscopy and virtual reality enabled us to reconstruct the capillary networks of twelve murine glomeruli and trace spatial graph models. We used circuit analysis to represent these as Reeb graphs, the fundamental theorem of calculus to describe a mean trajectory and its curvature, and principal components analysis to reveal lateral and transverse symmetry. Separately, we built a non-spatial random graph growth model based on two mechanisms, angiogenesis and intussusception, which provided evidence that both contribute to development. We then introduced several topological measures of lobularity and found, surprisingly, that empirical glomeruli tend to be less lobular than those generated by our model. Ongoing work focuses on simulation-based attack tolerance and the development of a spatial growth model.

Recently, the Nobel Prize winning machine learning (ML) model, AlphaFold, expanded its protein structure prediction capabilities from monomers to multimers with AlphaFold3. Here, we investigate this expanded utility in the novel context of mechanistic models of cell signaling. These models describe cell signaling events, such as binding, amongst a network of molecules, mostly proteins, and have been applied to answer both clinical and fundamental biology questions. For example, cell signaling models have been used to propose improvements to CAR-T cell therapies and to elucidate cellular ‘decision making’. Use of these models is oftentimes limited by a sparsity of data for parameterization. Thus, in this work, we introduce a Bayesian framework that incorporates information about protein structure to guide parameter inference for mechanistic models. Rather than searching all plausible parameter values, we can refine our search by considering information that is specific to the proteins involved in the signaling event. Excitingly, we find augmenting mechanistic models of signaling to be uniquely compatible with established ML models. We test our approach on two signaling models. In both cases, our approach improves parameter estimates—however, these improvements do not significantly change prediction performance on test data for either model. We find that this is due to a lack of sensitivity between the informed parameters and the test outputs. In contrast, when we examine an output that is sensitive to the changed parameters, we see a clear change in the predicted dynamics. We note that our proposed approach is limited to parameters of reversible, bimolecular binding reactions. Yet, excitingly, mechanistic models of cell signaling are often comprised of such reactions, ensuring the relatively wide applicability of our inference approach in this context.

Most antibiotics kill bacteria by disrupting cell wall formation during mitosis. Bacterial persisters are individuals within a population that avoid this fate by not replicating. We use a parabolic PDE to model the phenotypic switch between normal, active bacteria and persisters along with a nonlocal birth-jump process that captures epigenetic inheritance. In addition, we relate bacterial population development to resource dynamics in order to depict a more realistic bacterial growth limit. Mathematically, the model consists of a non-local PDE coupled to an ODE. We prove the well-posedness of the model using semi-group theory and the Banach fixed point theorem. We then examine the evolutionarily stable strategies of persister cells by conducting a global invasion analysis with an appropriately chosen Lyapunov functional.

The cornea represents the outermost transparent layer of the eye and is structured in multiple layers. The corneal epithelium, which forms the exterior surface of the cornea, is distinct from many other epithelial tissues in that it consists of 5 to 7 layers, rather than a single layer. This stratification process involves the upward movement of cells from the basal layer to the upper layers, a mechanism known as cell delamination. Additionally, the integrity of the corneal epithelium is maintained through the migration of new basal cells from the periphery toward the centre of the cornea. Despite its crucial role in maintaining corneal function, the regulatory mechanisms governing this process, as well as how it adapts to cell loss during wound healing, remain poorly understood. Our research aims to explore the regulation of corneal cell behavior through the use of a Voronoi cell-based model, which links local cellular interactions to the emergent dynamics observed in the stratified epithelium.

In this talk, I will illustrate how to model the dynamics of the fluorophores used in single-molecule localization microscopy (SMLM) as a hidden Markov chain with discretized emissions. I will generalize the proposed models in literature into a simple framework model. With the 3-state model as a particular example of our general formulation, I will show the process to obtain the transmission matrix by constructing a system of linear inhomogeneous transport partial differential equations (PDEs), which is solved by repeated Laplace-Inverse Laplace transforms. To demonstrate the usefulness of the transmission matrix, we designed two simple algorithms to solve the inference problem of the transition rates. In conclusion, the general formation is widely applicable to various techniques in SMLM, representative of the SMLM camera and adaptable to solve other active research problems such as molecule counting problems.

Proteins are macromolecules in the cell performing most of the metabolic processes. The protein is made up of a linear chain of amino-acids (primary structure) synthesized, through transcription and translation of the corresponding gene/DNA sequence inside the cell. The functional protein is a three-dimensional structure that is formed due to spontaneous or assisted folding of the linear chain decided by the physicochemical forces exerted due to the size, charge and chemical nature of the amino acids. The 3D structure essentially determines the function of the protein - known as the 'Structure-Function paradigm' in molecular biophysics. We have modelled the 3-dimensional structure of proteins using the network/graph theory, where the amino acids are the nodes, and links are the physicochemical forces that hold any two amino acids together. I will show how the network approach can clearly explain the large functional differences in proteins and their mutants, having insignificant structural variations, not easily identifiable using standard structural biology methods, and thereby questioning the universality of the 'Structure-Function paradigm'.

Cell signaling is governed by complex networks of biochemical interactions. These networks are critical for a wide range of cellular functions, including detecting environmental changes and cellular motility. Modeling these processes with reaction-diffusion equations (RDEs) requires prior knowledge of protein-protein interactions for constructing the underlying network. The complex nature of signaling pathways means many relevant interactions may be unknown. To address this challenge, we developed a deep learning-based method to infer reaction networks from data. By integrating a physics-informed neural network (PINN) with a neural network for symbolic regression, this method learns interpretable RDE models from spatiotemporal data, effectively learning the biochemical reactions driving dynamics. To develop and validate our approach, we applied it to data generated from a model of cell polarity establishment. This approach has the potential to overcome limitations from incomplete knowledge of protein-protein interactions, serving as a powerful tool for uncovering how cells regulate complex behaviors.

The stability and polarity of the microtubule cytoskeleton is required for long-range, sustained transport within neuronal cellsl. In particular, the healthy microtubule cytoskeleton is comprised of tubulin protein and is stable with a particular orientation. However, when injured, these microtubules are dynamic, rearrange their orientation, and the appearance of microtubules is upregulated. It is unknown what mechanisms are involved in this balance between dynamic rearrangement and sustained function. Using a stochastic mathematical model that incorporates experimental data, we seek to understand how nucleation can impact microtubule dynamics in dendrites of fruit fly neurons. In the stochastic model, we assume two mechanisms limit microtubule growth: limited tubulin availability and the dependence of shrinking events on microtubule length. To better understand our stochastic model, we develop a partial differential equation (PDE) model that describes microtubule growth and nucleation dynamics, and we compare analytical results to results from the complex stochastic model. Insights from these models can then be used to understand what mechanisms are used organize into polarized structures in neurons, and how microtubule dynamics, like nucleation, may impact cargo localization post-injury.

Tissue morphogenesis involves the self-organized creation of patterns and shapes. In many cases details of underlying mechanisms are elusive, yet an increasing amount of experimental data suggests that chemical morphogens and mechanical processes are strongly coupled. Here, we develop and simulate a minimal model for the emergence of asymmetry in aggregates of the Hydra polyp based on mechanochemical coupling of surface stiffness and a morphogen concentration. We contrast this model with the classical morphogen patterning mechanisms based on Turing type reaction diffusion systems. In analogy to this classical mechnism, we carry out the stability analysis of the lower dimensional toy model and identify minimal conditions for symmetry breaking. Our results suggest that mechanochemical pattern formation underlies symmetry breaking in Hydra.

Ribosomal stalling during translation of mRNA can result for example from oxidative conditions surrounding the site of translation. It impacts the cellular protein production machinery and therefore decrease cell proliferation. Accordingly, the rate of protein synthesis (R) can be considered as a hallmark of ribosomal stalling. In vivo experiments can determine protein content of a cell and differences in ribosomal density for different stalling scenarios. We employ the Ribosome Flow Model (RFM) coupled with Monte Carlo simulations to quantitatively establish the implications of three stalling patterns motivated by biological processes: We consider (1) the overall frequency of stalling sides as defined by harmful mRNA modification, (2) the degree of the reduction of the translocation rate λ reflecting the severity of mRNA transcriptional impairments, as well as (3) the effect of clustering, chain and gap impairments, as well as cluster locality of these anomalies. Each of these stalling patterns impacts protein synthesis rate and ribosomal density differently. We show how quantitative prediction of the impact of each and combinations of these patterns can be used as to study and predict mRNA stalling. Major findings of our analysis are, that for a given severity of mRNA damage, the equilibrium rate of protein synthesis R* does not depend on impairment locality, and is not related to the ribosomal density. In contrast, ribosomal density is strongly dependent on the locality of impairment clustering.

How pluripotent cells give rise to progressively more specialised cells over multiple cell divisions, known as cell fate, remains one of the mysteries of systems biology. During development, it is of the utmost importance that cells uphold certain division regimes for an organism to survive and thrive. Beyond development, cell fate perturbations can result in cancer and other pathological conditions. The theoretical and mathematical biology community has been making contributions to our understanding of cell fate including by quantifying Waddington’s seminal landscape using dynamical systems, performing statistical trajectory inference on single-cell sequencing data, or considering geometric and algebraic approaches to cell fate. In this talk, I will present a random dynamical systems interpretation of cell fate. This approach is, arguably, a generalisation of existing models of cell fate that may be able to provide new perspectives into cell fate.

Advances in cell type and gene expression mapping have significantly enhanced our understanding of the human body. However, comprehending interactions at cellular and tissue levels is equally critical for unveiling mechanisms underlying health and aging. This project leverages data from the Human BioMolecular Atlas Program (HuBMAP) and Human Cell Atlas (HCA) to explore cellular functionality within functional tissue units (FTUs) of the skin epidermis. Using HuBMAP and HCA single-cell RNA sequencing (scRNA-seq) and transcriptomics data, we identify key cell types acting as chemical secretors (ligands) and receivers (receptors) in healthy and diseased tissue. We employ PhysiCell, an agent-based modeling platform, to construct a 3D computational cellular environment. The workflow involves preprocessing the transcriptomics data into a machine-readable format and generating chemical communication graphs that capture th ofe dynamic interplay signaling molecules between secretor and receiver cells. By combining biological data with multiscale ABM, we aim to visually and quantitatively model the chemical interactions within the epidermal FTUs of human skin tissue. The overarching goal is to develop a mathematical model elucidating how hydration-mediated cellular communication impacts tissue homeostasis and delays aging processes. This research has the potential to provide new insights into the mechanisms of skin aging and inform strategies for promoting tissue health through hydration management.

In eukaryotic chemotaxis, cells extend and retract transient actin-driven protrusions at their membrane. These protrusions facilitate both the detection of external chemical gradients and directional movement via the formation of focal adhesions with the extracellular matrix. While extensive experimental work has characterised how protrusive activity varies with a range of environmental parameters, the mechanistic principles governing these relationships remain poorly understood. Here, we model the extension of actin-based protrusions in chemotaxis mathematically as an optimisation problem, wherein cells must balance the detection of external chemical gradients with the energetic cost of protrusion formation. The model highlights energetic efficiency in movement as a major predictor of phenotypic variation amongst motile cell populations, successfully reproducing experimentally observed but previously non-understood patterns of protrusive activity across a range of biological systems. Additionally, we leverage the model to generate novel predictions regarding cellular responses to environmental perturbations, providing testable hypotheses for future experimental work.

Cell signaling is governed by complex networks of biochemical interactions. These networks are critical for a wide range of cellular functions, including detecting environmental changes and cellular motility. Modeling these processes with reaction-diffusion equations (RDEs) requires prior knowledge of protein-protein interactions for constructing the underlying network. The complex nature of signaling pathways means many relevant interactions may be unknown. To address this challenge, we developed a deep learning-based method to infer reaction networks from data. By integrating a physics-informed neural network (PINN) with a neural network for symbolic regression, this method learns interpretable RDE models from spatiotemporal data, effectively learning the biochemical reactions driving dynamics. To develop and validate our approach, we applied it to data generated from a model of cell polarity establishment. This approach has the potential to overcome limitations from incomplete knowledge of protein-protein interactions, serving as a powerful tool for uncovering how cells regulate complex behaviors.

The stability and polarity of the microtubule cytoskeleton is required for long-range, sustained transport within neuronal cellsl. In particular, the healthy microtubule cytoskeleton is comprised of tubulin protein and is stable with a particular orientation. However, when injured, these microtubules are dynamic, rearrange their orientation, and the appearance of microtubules is upregulated. It is unknown what mechanisms are involved in this balance between dynamic rearrangement and sustained function. Using a stochastic mathematical model that incorporates experimental data, we seek to understand how nucleation can impact microtubule dynamics in dendrites of fruit fly neurons. In the stochastic model, we assume two mechanisms limit microtubule growth: limited tubulin availability and the dependence of shrinking events on microtubule length. To better understand our stochastic model, we develop a partial differential equation (PDE) model that describes microtubule growth and nucleation dynamics, and we compare analytical results to results from the complex stochastic model. Insights from these models can then be used to understand what mechanisms are used organize into polarized structures in neurons, and how microtubule dynamics, like nucleation, may impact cargo localization post-injury.

Tissue morphogenesis involves the self-organized creation of patterns and shapes. In many cases details of underlying mechanisms are elusive, yet an increasing amount of experimental data suggests that chemical morphogens and mechanical processes are strongly coupled. Here, we develop and simulate a minimal model for the emergence of asymmetry in aggregates of the Hydra polyp based on mechanochemical coupling of surface stiffness and a morphogen concentration. We contrast this model with the classical morphogen patterning mechanisms based on Turing type reaction diffusion systems. In analogy to this classical mechnism, we carry out the stability analysis of the lower dimensional toy model and identify minimal conditions for symmetry breaking. Our results suggest that mechanochemical pattern formation underlies symmetry breaking in Hydra.

Ribosomal stalling during translation of mRNA can result for example from oxidative conditions surrounding the site of translation. It impacts the cellular protein production machinery and therefore decrease cell proliferation. Accordingly, the rate of protein synthesis (R) can be considered as a hallmark of ribosomal stalling. In vivo experiments can determine protein content of a cell and differences in ribosomal density for different stalling scenarios. We employ the Ribosome Flow Model (RFM) coupled with Monte Carlo simulations to quantitatively establish the implications of three stalling patterns motivated by biological processes: We consider (1) the overall frequency of stalling sides as defined by harmful mRNA modification, (2) the degree of the reduction of the translocation rate λ reflecting the severity of mRNA transcriptional impairments, as well as (3) the effect of clustering, chain and gap impairments, as well as cluster locality of these anomalies. Each of these stalling patterns impacts protein synthesis rate and ribosomal density differently. We show how quantitative prediction of the impact of each and combinations of these patterns can be used as to study and predict mRNA stalling. Major findings of our analysis are, that for a given severity of mRNA damage, the equilibrium rate of protein synthesis R* does not depend on impairment locality, and is not related to the ribosomal density. In contrast, ribosomal density is strongly dependent on the locality of impairment clustering.

How pluripotent cells give rise to progressively more specialised cells over multiple cell divisions, known as cell fate, remains one of the mysteries of systems biology. During development, it is of the utmost importance that cells uphold certain division regimes for an organism to survive and thrive. Beyond development, cell fate perturbations can result in cancer and other pathological conditions. The theoretical and mathematical biology community has been making contributions to our understanding of cell fate including by quantifying Waddington’s seminal landscape using dynamical systems, performing statistical trajectory inference on single-cell sequencing data, or considering geometric and algebraic approaches to cell fate. In this talk, I will present a random dynamical systems interpretation of cell fate. This approach is, arguably, a generalisation of existing models of cell fate that may be able to provide new perspectives into cell fate.

Advances in cell type and gene expression mapping have significantly enhanced our understanding of the human body. However, comprehending interactions at cellular and tissue levels is equally critical for unveiling mechanisms underlying health and aging. This project leverages data from the Human BioMolecular Atlas Program (HuBMAP) and Human Cell Atlas (HCA) to explore cellular functionality within functional tissue units (FTUs) of the skin epidermis. Using HuBMAP and HCA single-cell RNA sequencing (scRNA-seq) and transcriptomics data, we identify key cell types acting as chemical secretors (ligands) and receivers (receptors) in healthy and diseased tissue. We employ PhysiCell, an agent-based modeling platform, to construct a 3D computational cellular environment. The workflow involves preprocessing the transcriptomics data into a machine-readable format and generating chemical communication graphs that capture th ofe dynamic interplay signaling molecules between secretor and receiver cells. By combining biological data with multiscale ABM, we aim to visually and quantitatively model the chemical interactions within the epidermal FTUs of human skin tissue. The overarching goal is to develop a mathematical model elucidating how hydration-mediated cellular communication impacts tissue homeostasis and delays aging processes. This research has the potential to provide new insights into the mechanisms of skin aging and inform strategies for promoting tissue health through hydration management.

In eukaryotic chemotaxis, cells extend and retract transient actin-driven protrusions at their membrane. These protrusions facilitate both the detection of external chemical gradients and directional movement via the formation of focal adhesions with the extracellular matrix. While extensive experimental work has characterised how protrusive activity varies with a range of environmental parameters, the mechanistic principles governing these relationships remain poorly understood. Here, we model the extension of actin-based protrusions in chemotaxis mathematically as an optimisation problem, wherein cells must balance the detection of external chemical gradients with the energetic cost of protrusion formation. The model highlights energetic efficiency in movement as a major predictor of phenotypic variation amongst motile cell populations, successfully reproducing experimentally observed but previously non-understood patterns of protrusive activity across a range of biological systems. Additionally, we leverage the model to generate novel predictions regarding cellular responses to environmental perturbations, providing testable hypotheses for future experimental work.

Forest harvesting in coastal British Columbia (BC) has altered watersheds, impacting salmon habitat by increasing sediment, reducing riparian cover, and altering hydrology. These changes can affect the survival and growth of salmon through mechanisms like reduced egg-to-fry survival, increased stream scour, increased thermal stress, and loss of stream complexity which is essential for salmon rearing. Despite numerous localized studies, no comprehensive analysis has examined the population-level effects of forestry on salmon across BC. After assembling forest harvest histories for 1,746 salmon-bearing watersheds (1883-2022) and salmon abundance data (1950-2022), we used stock-recruit models (Ricker and Beverton-Holt) in a hierarchical Bayesian framework to test the effects of forestry metrics (Equivalent Clearcut Area, Cumulative Percent Disturbed) on chum and pink salmon productivity. Our results show a strong negative effect of forestry on chum productivity (e.g. 25% equivalent clearcut area reduces productivity by more than 20%) and a negligible effect on pink salmon. This highlights forestry’s significant role in the decline of chum salmon populations over recent decades.

Due to their diversity and abundance, insects play essential ecological roles, including crop pollination, nutrient cycling, and serving as a food source for other species. However, climate change is predicted to heavily impact insect populations, with some expected to decline by up to 18% globally by the end of the 2020s, raising concerns about the future health of the bioeconomy. To investigate these impacts, we conducted a temperature study on Aster leafhoppers (Hemiptera: Cicadellidae: Macrosteles quadrilineatus). Using five temperature conditions, we collected physiological data over a month to assess the impact of temperature on survival, maturation, and egg production. We then developed an age-structured population model to explore how environmental temperature influences insect fitness and mortality rates. Our model is parameterized with experimental data across various climate change scenarios, providing insights into the effects of rising temperatures on insect survival and population dynamics. This work highlights the cascading effects of climate change on ecological networks and emphasizes the importance of understanding insect responses to environmental stressors.

Recent progress in genome assembly techniques has led to an explosion in chromosome-length genome sequences. These unfragmented assemblies have enabled biologists to study molecular evolution at unprecedented scales, providing insight into the evolution of genome architecture. Microsynteny, the conservation of gene order, has proven to be a key concept in our understanding of genome evolution. However, it remains unclear when microsynteny occurs due to random chance or selection. Here, we develop a mathematical model to discriminate between these two cases. Our model describes the dynamics of synteny block size distributions in the absence of selection or other biases. By fitting this null model to data from a comparative analysis of mammalian genomes, we identify synteny blocks larger than expected in the absence of selection. This approach allows us to rigorously determine which sets of genes are likely to have selection on their ordering in a lineage-specific manner. Our model presents a powerful tool for uncovering functional relationships between genes based on their ordering and for understanding the evolution of gene co-regulation.

Empirical evidence shows that evolution may take place during species' range expansion. Indeed, dispersal ability tends to be selected for at the leading edge of invasions, ultimately increasing a species' spreading speed. However, for organisms across many different taxa, higher dispersal comes at the cost of fitness, producing evolutionary trade-offs at the leading edge. Using reaction-diffusion equations and adaptive dynamics, we provide new insights on how such evolutionary processes take place. We show how evolution may drive phenotypes at the leading edge to maximize the asymptotic spreading speed, and conditions under which phenotypic plasticity in dispersal is selected for under different dispersal-reproduction trade-off scenarios. We provide some possible future research directions and other systems where the framework can be applied.

Cyanobacterial blooms (CBs) pose significant global challenges due to their harmful toxins and socio-economic impacts, with nutrient availability playing a key role in their growth, as described by ecological stoichiometry (ES). However, real-world ecosystems exhibit spatial heterogeneity, limiting the applicability of simpler, spatially uniform models. To address this, we develop a spatially explicit partial differential equation model based on ES to study cyanobacteria in the epilimnion of freshwater systems. We establish the well-posedness of the model and perform a stability analysis, showing that it admits two linearly stable steady states, leading to either extinction or saturation. We use the finite elements method to numerically solve our system on a real lake domain derived from Geographic Information System (GIS) data and realistic wind conditions extrapolated from ERA5-Land. Our numerical results highlight the importance of lake shape and size in CB monitoring, while global sensitivity analysis using Sobol Indices identifies light attenuation and intensity as primary drivers of bloom variation, with water movement influencing early bloom stages and nutrient input becoming critical over time. This model supports continuous water-quality monitoring, informing agricultural, recreational, economic, and public health strategies for mitigating CBs.

Whether interspecific competition is a major contributing factor in setting species' range limits has been debated for a long time. Theoretical studies using evolutionary models have proposed that the interaction between interspecific competition and disruptive gene flow along an environmental gradient can halt range expansion of ecologically related species where they meet. However, the stability of such range limits has not been well addressed. In this talk, I present our work on investigating the stability of competitively formed range limits using a deterministic model of adaptive range evolution. We show that the range limits are unlikely to be evolutionarily stable if the environmental optima for fitness-related traits vary linearly in space. However, we demonstrate that environmental nonlinearities such as “knees” and “wiggles”, wherein an isolated sharp change or a step-like change occurs in the steepness of a trait optimum, can strongly stabilize the range limits. We show that the stability of the range limits established at such nonlinearities is robust against moderate environmental disturbances. Although strong climatic changes can still destabilize the range limits, such destabilization depends on how the relative dominance of the competing species changes across the environmental nonlinearity. Therefore, our results highlight the importance of measuring the competitive ability of species when predicting their response to climate change.

Tumours are complex ecosystems where diverse cancer cell subpopulations interact with each other and with non-cancer cells around them. Evolutionary game theory (EGT) has established itself as a powerful mathematical framework to study the implications of such ecological interactions, demonstrating an important role in shaping oncogenesis and treatment response. However, much of this work has been theoretical using parameters that are only loosely grounded in biological data. To move towards quantitative and predictive models of tumour ecology it is crucial to develop theoretical and experimental methodology to empirically calibrate and validate EGT models. We present an in silico study to optimize the 'Game Assay' for measuring ecological interactions between cancer cell populations in vitro. This assay, originally developed by Kaznatcheev et al (2017), involves co-culturing populations at different ratios, monitoring growth rates via time-lapse microscopy, and inferring frequency-dependent interactions. We begin by characterizing the accuracy and precision of this assay in a simulation study in which we use the replicator equation as the “ground truth”. Our simulations reveal potential biases in estimating fitness differences and interaction parameters, highlighting the need for careful experimental design. We provide guidelines for optimizing seeding ratios, number of replicates, and frequency of measurements, and present a new analysis techniques to improve the accuracy and precision of interaction measurements. Finally, we apply our optimized protocol to quantify interactions between 4 drug-sensitive and resistant lung cancer cell lines, revealing diverse ecological dynamics. This work demonstrates the power of integrating mathematical modeling with experimental approaches to develop robust empirical protocols and gain a quantitative understanding of tumour ecology.

Genetically engineered mosquitoes (GEMs) offer a promising malaria control strategy, yet their ecological interactions, dispersal, and long-term effects remain uncertain. Accurate modeling is essential to optimize GEM release strategies and assess their effectiveness in natural ecosystems. This study presents a high-performance, exascale agent-based model (ABM) simulating gene drive dynamics in wild mosquito populations. Incorporating mosquito population dynamics, spatial ecology, and genotype inheritance, the model provides insights into optimizing release timing, locations, and dispersal strategies. Our findings indicate that under optimal dispersal conditions, GEMs can achieve a 95% prevalence in wild populations within 112 days. Furthermore, our findings indicate that strategically coordinating GEM releases across multiple sites does not significantly impact gene drive establishment on the island. By capturing mosquito behaviors and movement in heterogeneous environments, this ABM serves as a powerful tool for evaluating GEM interventions, supporting evidence-based malaria control strategies, and enhancing ecological understanding of gene drive propagation..

Monitoring cyanobacterial blooms promptly and accurately is crucial for public health management and understanding aquatic ecosystem dynamics. Remote sensing, particularly satellite observations, offers a viable approach for continuous monitoring. This study utilizes multispectral images from the Sentinel-2 satellite constellation in conjunction with ERA5-Land data to facilitate broad-scale data collection. We proposed a simple deep convolutional neural network (CNN) architecture to analyze cyanobacteria (CB) concentration dynamics in Pigeon Lake, Canada, over a five-year period. Utilizing the Local Getis-Ord statistic, we identified and analyzed trends in hot and cold spots under the null hypothesis of random distribution. We observed changes in the distribution and median CB concentration in hot spots over time. Additionally, we trained a Kolmogorov-Arnold Network (KAN) to classify segments of the lake shoreline into hot and non-hot spots using the Dynamic World dataset within a 500-meter radius of the lake.

Microalgae biofilms have been demonstrated to recover nutrients from wastewater and serve as biomass feedstock for bioproducts. However, there is a need to develop a platform to quantitatively describe microalgae biofilm production, which can provide guidance and insights for improving biomass areal productivity and nutrient uptake efficiency. In this talk, I will introduce a unified experimental and theoretical framework to investigate algae biofilm growth on a rotating algae biofilm reactor (RABR). Experimental laboratory setups are used to conduct controlled experiments on testing environmental and operational factors for RABRs. We propose a differential–integral equation‐based mathematical model for microalgae biofilm cultivation guided by laboratory experimental findings. The predictive mathematical model development is coordinated with laboratory experiments of biofilm areal productivity associated with ammonia and inorganic phosphorus uptake by RABRs. The unified experimental and theoretical tool is used to investigate the effects of RABR rotating velocity, duty cycle (DC), and light intensity on algae biofilm growth, areal productivity, nutrient uptake efficiency, and energy efficiency in wastewater treatment.

Mosquitoes are key vectors of several infectious diseases affecting humans and animals. In North America, Culex mosquitoes are primary vectors of West Nile virus, St. Louis encephalitis, and Japanese encephalitis, as well as viral diseases in birds and horses. The Culex mosquito life cycle consists of four stages: eggs, larvae, pupae, and adults, each with unique development and mortality rates. Only active (non-diapausing) adults can reproduce, and environmental factors such as temperature, photoperiod, and rainfall influence population dynamics and stage-specific abundances. We develop a data-driven, stage-structured model that incorporates experimental data to describe how key climate variables regulate life history parameters. Specifically, egg laying rates depend on temperature, while maturation and survival rates are influenced by both temperature and rainfall. Mortality is temperature-dependent, and diapause induction and reactivation rates in adults are driven by temperature and photoperiod. Unlike many previous models that focus on tropical mosquitoes, our study explicitly includes diapause, a dormancy period in adult Culex mosquitoes essential for accurate modelling of temperate mosquito populations. Our results show that mosquito populations peak during summer months when temperatures exceed 10°C. Seasonal fluctuations in abundance highlight the need for adaptive vector control strategies. Since control measures often target specific life stages, such as larvicides for larvae or insecticides for adults, our findings suggest that optimal intervention strategies should vary by season to effectively reduce mosquito populations and disease risk.

Many biological systems regulate phenotypic heterogeneity as a fitness-maximising strategy in uncertain and dynamic environments. Analysis of such strategies is typically confined both to a discrete set of environmental conditions, and to a discrete (often binary) set of phenotypes specialised to each condition. In this talk, we extend on both fronts to encapsulate both a discrete and continuous spectrum of phenotypes arising in response to two broad classes of environmental fluctuations that drive changes in the phenotype-dependent growth rates. We present a series of analytical and semi-analytical results that reveal regimes in which both discrete and continuous phenotypic heterogeneity is evolutionary advantageous.

In this study, we examine how explicit spatial memory influences different mathematical models in various ecological dispersal contexts. Specifically, we analyze three memory-based dispersal strategies: (1) gradient-based movement, where individuals respond to environmental gradients; (2) environment matching, which promotes a uniform population distribution; and (3) location-based movement, where decisions are based solely on local suitability. These strategies correspond to diffusion-advection, Fickian diffusion, and Fokker-Planck diffusion models, respectively. Additionally, we explore steady-state problems for each strategy to highlight the differences between models incorporating temporal memory and those incorporating spatiotemporal memory.

Whether interspecific competition is a major contributing factor in setting species' range limits has been debated for a long time. Theoretical studies using evolutionary models have proposed that the interaction between interspecific competition and disruptive gene flow along an environmental gradient can halt range expansion of ecologically related species where they meet. However, the stability of such range limits has not been well addressed. In this talk, I present our work on investigating the stability of competitively formed range limits using a deterministic model of adaptive range evolution. We show that the range limits are unlikely to be evolutionarily stable if the environmental optima for fitness-related traits vary linearly in space. However, we demonstrate that environmental nonlinearities such as “knees” and “wiggles”, wherein an isolated sharp change or a step-like change occurs in the steepness of a trait optimum, can strongly stabilize the range limits. We show that the stability of the range limits established at such nonlinearities is robust against moderate environmental disturbances. Although strong climatic changes can still destabilize the range limits, such destabilization depends on how the relative dominance of the competing species changes across the environmental nonlinearity. Therefore, our results highlight the importance of measuring the competitive ability of species when predicting their response to climate change.

Tumours are complex ecosystems where diverse cancer cell subpopulations interact with each other and with non-cancer cells around them. Evolutionary game theory (EGT) has established itself as a powerful mathematical framework to study the implications of such ecological interactions, demonstrating an important role in shaping oncogenesis and treatment response. However, much of this work has been theoretical using parameters that are only loosely grounded in biological data. To move towards quantitative and predictive models of tumour ecology it is crucial to develop theoretical and experimental methodology to empirically calibrate and validate EGT models. We present an in silico study to optimize the 'Game Assay' for measuring ecological interactions between cancer cell populations in vitro. This assay, originally developed by Kaznatcheev et al (2017), involves co-culturing populations at different ratios, monitoring growth rates via time-lapse microscopy, and inferring frequency-dependent interactions. We begin by characterizing the accuracy and precision of this assay in a simulation study in which we use the replicator equation as the “ground truth”. Our simulations reveal potential biases in estimating fitness differences and interaction parameters, highlighting the need for careful experimental design. We provide guidelines for optimizing seeding ratios, number of replicates, and frequency of measurements, and present a new analysis techniques to improve the accuracy and precision of interaction measurements. Finally, we apply our optimized protocol to quantify interactions between 4 drug-sensitive and resistant lung cancer cell lines, revealing diverse ecological dynamics. This work demonstrates the power of integrating mathematical modeling with experimental approaches to develop robust empirical protocols and gain a quantitative understanding of tumour ecology.

Genetically engineered mosquitoes (GEMs) offer a promising malaria control strategy, yet their ecological interactions, dispersal, and long-term effects remain uncertain. Accurate modeling is essential to optimize GEM release strategies and assess their effectiveness in natural ecosystems. This study presents a high-performance, exascale agent-based model (ABM) simulating gene drive dynamics in wild mosquito populations. Incorporating mosquito population dynamics, spatial ecology, and genotype inheritance, the model provides insights into optimizing release timing, locations, and dispersal strategies. Our findings indicate that under optimal dispersal conditions, GEMs can achieve a 95% prevalence in wild populations within 112 days. Furthermore, our findings indicate that strategically coordinating GEM releases across multiple sites does not significantly impact gene drive establishment on the island. By capturing mosquito behaviors and movement in heterogeneous environments, this ABM serves as a powerful tool for evaluating GEM interventions, supporting evidence-based malaria control strategies, and enhancing ecological understanding of gene drive propagation..

Monitoring cyanobacterial blooms promptly and accurately is crucial for public health management and understanding aquatic ecosystem dynamics. Remote sensing, particularly satellite observations, offers a viable approach for continuous monitoring. This study utilizes multispectral images from the Sentinel-2 satellite constellation in conjunction with ERA5-Land data to facilitate broad-scale data collection. We proposed a simple deep convolutional neural network (CNN) architecture to analyze cyanobacteria (CB) concentration dynamics in Pigeon Lake, Canada, over a five-year period. Utilizing the Local Getis-Ord statistic, we identified and analyzed trends in hot and cold spots under the null hypothesis of random distribution. We observed changes in the distribution and median CB concentration in hot spots over time. Additionally, we trained a Kolmogorov-Arnold Network (KAN) to classify segments of the lake shoreline into hot and non-hot spots using the Dynamic World dataset within a 500-meter radius of the lake.

Microalgae biofilms have been demonstrated to recover nutrients from wastewater and serve as biomass feedstock for bioproducts. However, there is a need to develop a platform to quantitatively describe microalgae biofilm production, which can provide guidance and insights for improving biomass areal productivity and nutrient uptake efficiency. In this talk, I will introduce a unified experimental and theoretical framework to investigate algae biofilm growth on a rotating algae biofilm reactor (RABR). Experimental laboratory setups are used to conduct controlled experiments on testing environmental and operational factors for RABRs. We propose a differential–integral equation‐based mathematical model for microalgae biofilm cultivation guided by laboratory experimental findings. The predictive mathematical model development is coordinated with laboratory experiments of biofilm areal productivity associated with ammonia and inorganic phosphorus uptake by RABRs. The unified experimental and theoretical tool is used to investigate the effects of RABR rotating velocity, duty cycle (DC), and light intensity on algae biofilm growth, areal productivity, nutrient uptake efficiency, and energy efficiency in wastewater treatment.

Mosquitoes are key vectors of several infectious diseases affecting humans and animals. In North America, Culex mosquitoes are primary vectors of West Nile virus, St. Louis encephalitis, and Japanese encephalitis, as well as viral diseases in birds and horses. The Culex mosquito life cycle consists of four stages: eggs, larvae, pupae, and adults, each with unique development and mortality rates. Only active (non-diapausing) adults can reproduce, and environmental factors such as temperature, photoperiod, and rainfall influence population dynamics and stage-specific abundances. We develop a data-driven, stage-structured model that incorporates experimental data to describe how key climate variables regulate life history parameters. Specifically, egg laying rates depend on temperature, while maturation and survival rates are influenced by both temperature and rainfall. Mortality is temperature-dependent, and diapause induction and reactivation rates in adults are driven by temperature and photoperiod. Unlike many previous models that focus on tropical mosquitoes, our study explicitly includes diapause, a dormancy period in adult Culex mosquitoes essential for accurate modelling of temperate mosquito populations. Our results show that mosquito populations peak during summer months when temperatures exceed 10°C. Seasonal fluctuations in abundance highlight the need for adaptive vector control strategies. Since control measures often target specific life stages, such as larvicides for larvae or insecticides for adults, our findings suggest that optimal intervention strategies should vary by season to effectively reduce mosquito populations and disease risk.

Many biological systems regulate phenotypic heterogeneity as a fitness-maximising strategy in uncertain and dynamic environments. Analysis of such strategies is typically confined both to a discrete set of environmental conditions, and to a discrete (often binary) set of phenotypes specialised to each condition. In this talk, we extend on both fronts to encapsulate both a discrete and continuous spectrum of phenotypes arising in response to two broad classes of environmental fluctuations that drive changes in the phenotype-dependent growth rates. We present a series of analytical and semi-analytical results that reveal regimes in which both discrete and continuous phenotypic heterogeneity is evolutionary advantageous.

In this study, we examine how explicit spatial memory influences different mathematical models in various ecological dispersal contexts. Specifically, we analyze three memory-based dispersal strategies: (1) gradient-based movement, where individuals respond to environmental gradients; (2) environment matching, which promotes a uniform population distribution; and (3) location-based movement, where decisions are based solely on local suitability. These strategies correspond to diffusion-advection, Fickian diffusion, and Fokker-Planck diffusion models, respectively. Additionally, we explore steady-state problems for each strategy to highlight the differences between models incorporating temporal memory and those incorporating spatiotemporal memory.

Defective interfering particles (DIPs) are mutated, replication-incompetent virions that can inhibit their corresponding standard virus (STV). Previous studies have shown the effectiveness of DIPs against various virus species, highlighting them as promising broad-spectrum antivirals. OP7, an influenza DIP with 37 nucleotide substitutions in its segment 7 (S7) vRNA, has been found to suppress STV replication more effectively than conventional DIPs. However, the effects of these mutations on the replication of OP7 and its mechanism of interference with the STV remained unclear. In this study, we investigated the infection dynamics during a coinfection of influenza STV and OP7 in cell culture. We monitored the dynamics of viral RNAs, assessed viral protein levels, and determined virus titers. With these experimental results, we developed a mathematical model to simulate the coinfection of STV and OP7. Subsequently, we used this model to explore various hypotheses about the impact of mutations on virus replication and to predict the suppression of STV by OP7 in passaging experiments. In vitro experiments show that S7-OP7 surpasses the levels of all STV genome segments. Model simulations suggest this is induced by a significantly increased rate of replication, attributed to mutations in S7-OP7 inducing a “superpromoter”. Additionally, simulations predicted a notable reduction in viral mRNA transcription for S7-OP7, which was later validated experimentally. Moreover, we deduce that the M1 protein derived from S7-OP7 mRNA is likely defective. Lastly, the model accurately predicts the spread of OP7 and the suppression of STV in infected cell cultures over multiple passages under various initial conditions. In summary, we developed a mathematical model that enables a thorough examination of STV and OP7 coinfection, improves our understanding of DIP interference mechanisms, and supports the development of antiviral therapies.

Chronic spontaneous urticaria (CSU) is a persistent skin disorder characterized by red, itchy eruptions of various shapes, known as wheals. These wheals appear and disappear daily, persisting for months or even decades, and severely impact patients' quality of life. The standard treatment for CSU primarily consists of second-generation H1 antihistamines, often administered at higher-than-usual doses. However, approximately 30% of patients remain symptomatic despite these conventional therapies. On the other hand, our previous mathematical modeling and clinical studies have identified five distinct types of wheal shapes through the development of clinical criteria for eruption geometry (EGe Criteria) and have shown how the characteristics of each wheal type are involved in the pathophysiology of CSU. These findings suggest that CSU may be classified into five medical subtypes based on wheal morphology. Thus, in this study, we explore the effectiveness of antihistamines based on wheal shape. We first evaluate the efficacy of antihistamines in silico using three key measures: wheal area, itching severity, and wheal expansion dynamics, across the five identified wheal types. Additionally, we validate some of our theoretical observations using clinical data from patients. By elucidating the relationships among the key networks involved in CSU pathophysiology, wheal morphology, and drug efficacy, we can enhance the development of more accurate diagnostic tools and treatment strategies in clinical settings.

Combination antiretroviral therapy (ART) can treat but not cure HIV, motivating the development of therapies that stimulate the immune system to control or eliminate infection. Two such immunotherapies- a TLR7 agonist and a therapeutic vaccine - were previously tested in SIV-infected rhesus macaques. Animals received ART alone or with concurrent single or combination immunotherapy, and viral rebound was monitored after treatment interruption. Many treated animals exhibited altered rebound kinetics, and a subset achieved either complete viral suppression or immune control after an initial rebound. However, the mechanisms driving these effects are unknown: do these therapies deplete the latent reservoir or enhance antiviral immunity, and do they act synergistically? To investigate the effects of immunotherapy, we built a mathematical model of viral dynamics incorporating latent cell reactivation and a generalized immune response. We confirmed the model could reproduce the range of rebound trajectories seen in the data, and examined whether parameters could be reliably estimated from the available data. Using nonlinear mixed-effects modeling, we quantified interindividual variability and identified significant differences in model parameters between treatment groups. Our results indicate that the vaccine alone reduces latent virus reactivation and enhances immune response avidity. The TLR7 agonist, when administered after late ART initiation, increases target cell availability and reduces the latent reservoir. We found that regardless of ART initiation, the two therapies act synergistically to further enhance immune response avidity. Immune avidity appeared to increase with later ART initiation, although whether this effect is specific to TLR7 treatment is unclear. Our model provides mechanistic insight into immunotherapeutic control of viral rebound and can be adapted to predict their impact in controlling HIV, guiding future therapeutic design and clinical trials.

The use of antibiotics during a disease outbreak presents a critical tradeoff between immediate treatment benefits to the individual and the long-term risk to the population. Typically, the extensive use of antibiotics has been thought to increase selective pressures, leading to resistance. This study explores scenarios where expanded antibiotic treatment can be advantageous for both individual and population health. We develop a mathematical framework to assess the impacts on outbreak dynamics of choosing to treat moderate infections not treated under current guidelines, focusing on cholera as a case study. We derive conditions under which treating moderate infections can sufficiently decrease transmission and reduce the total number of antibiotic doses administered. We identify two critical thresholds: the Outbreak Prevention Threshold (OPT), where expanded treatment reduces the reproductive number below one and halts transmission, and the Dose Utilization Threshold (DUT), where expanded treatment results in fewer total antibiotic doses used than under current guidelines. For cholera, we find that treating moderate infections can feasibly stop an outbreak when the untreated reproductive number is less than 1.42 and will result in fewer does used compared to current guidelines when the untreated reproductive number is less than 1.53. These findings demonstrate that conditions exist under which expanding treatment to include moderate infections can reduce disease spread and the selective pressure for antibiotic resistance. These findings extend to other pathogens and outbreak scenarios, suggesting potential targets for optimized treatment strategies that balance public health benefits and antibiotic stewardship.

Brazil has experienced recurrent dengue outbreaks, with over 18 million reported cases since 2000 and a record-breaking surge in 2024. Notably, there has been a demographic shift in disease burden, with an increasing proportion of severe cases and fatalities among the elderly. Current vaccination strategies, including the WHO-recommended Qdenga® rollout for children, may not effectively address this emerging risk. This study employs a mathematical modeling approach to evaluate age-specific and geographically targeted vaccination strategies. A Susceptible-Infected-Recovered (SIR)-based model, calibrated using historical dengue data from Brazil's Notifiable Diseases Information System (SINAN), incorporates spatial heterogeneity across 27 states and demographic factors such as prior exposure and birth rates. We assess the impact of different vaccination strategies by estimating averted infections, hospitalizations, fatalities, and years of life lost. Preliminary findings indicate significant variation in the force of infection across states and suggest that prioritizing vaccination for elderly populations may substantially reduce severe disease burden. This modeling framework provides a quantitative basis for optimizing vaccination policies, with potential applications to other arboviral diseases and endemic settings worldwide.

The emergence of infectious diseases remains a huge challenge to public health. Early detection of outbreaks using early warning signals (EWSs) offers an invaluable opportunity for effective preparedness and disease management. In this study, we seek to understand the structure of these signals using a mechanistic model that captures epidemic and social processes. We analyzed the simulated time series for change points and EWSs (autocorrelation and variance). All time series showed the expected delay in that the detected change point occurred significantly after the parameter passed the bifurcation point. These early warning signals exhibited a stronger response after the threshold for disease emergence had been exceeded. Assessing different disease progression and intervention models will help determine the most effective signals for use in public-health settings.

COVID-19, a global pandemic with severe health and economic repercussions, has prompted various approaches to mitigate its impact. We adapt an age-structured SEIVS model—incorporating waning immunity and partial protection—to explore counterfactual scenarios of non-pharmaceutical interventions (e.g., school/workplace closures, distancing) and selected vaccine policy changes in Canada. By altering contact patterns and compliance levels, we assess potential outcomes under stricter, earlier, or more relaxed mitigation measures, with a lesser emphasis on shifting vaccination rollouts. Findings indicate that timely, robust mitigation can substantially reduce severe disease and delay epidemic peaks, whereas delayed or minimal interventions lead to higher case burden. Integrating vaccine strategies into these scenarios further highlights the interplay between pharmaceutical and non-pharmaceutical measures, showcasing how modeling can inform dynamic policy-making for ongoing and future public health crises.

Advances in experimental techniques allow the collection of high-resolution spatio-temporal data that track individual motile entities over time and could be used to calibrate mathematical models of individual motility. However, experimental data is intrinsically discrete and noisy, and these characteristics complicate the effective calibration of models for individual motion. We consider individuals whose movement can be described by velocity-jump models in one spatial dimension, characterised by successive Markovian transitions between a network of n states, each with a specified velocity and a fixed rate of switching to every other state. We develop a Bayesian framework to calibrate these models to discrete and noisy data, which uses a likelihood consisting of approximations to the model solutions which we previously obtained. We apply the framework to recover the model parameters of simulated data, including the probabilities of switching to every other state. Moreover, we test the ability of the framework to select the most appropriate model to fit the data, including comparisons varying the number of states n.

During sleep periods, most mammals alternate multiple times between rapid-eye-movement (REM) sleep and non-REM (NREM) sleep. A common theory proposes that these transitions are governed by an ``hourglass-like'' homeostatic need to enter REM sleep that accumulates during the inter-REM interval and partially discharges during REM sleep. However, markers or mechanisms for REM homeostatic pressure remain undetermined. Recently, an analysis of sleep in mice demonstrated that the cumulative distribution function (CDF) of the amount of NREM sleep between REM bouts correlates with REM bout duration, suggesting that time in NREM sleep influences REM sleep need. Here, we build on those results and construct a predictive measure for the propensity to enter REM sleep as a function of time in NREM sleep since the previous REM episode. The REM propensity measure is precisely defined as the probability to enter REM sleep before the accumulation of an additional pre-specified amount of NREM sleep. Analyzing spontaneous sleep in mice, we find that, as NREM sleep accumulates between REM bouts, the REM propensity exhibits a peak value and then decays to zero with further NREM accumulation. We show that the REM propensity at REM onset predicts features of the subsequent REM bout under certain conditions. Specifically, during the light phase and for REM propensities occurring before the peak in propensity, the REM propensity at REM onset is correlated with REM bout duration, and with the probability of the occurrence of a short REM cycle called a sequential REM cycle. Further, we also find that proportionally more REM sleep occurs during sequential REM cycles, supporting a correlation between high values of our REM propensity measure and high REM sleep need. These results support the theory that a homeostatic need to enter REM sleep accrues during NREM sleep, but only for a limited range of NREM sleep accumulation. Time permitting, we will discuss current research directions.

CAR-T cell therapy, where patient T cells are genetically modified to target CD19-presenting B cells, has transformed the treatment landscape for several types of B-cell lymphomas. However, this therapy often triggers a strong inflammatory response, which can cause the tumor to temporarily swell in the days following CAR-T infusion — a phenomenon known as pseudo-progression. In this work, we present a dynamical model of CAR-T cell therapy that explicitly incorporates the role of pro-inflammatory cytokines in shaping treatment outcomes. Our model reproduces a wide spectrum of clinical trajectories, including complete remission, treatment failure, and transient pseudo-progression. Importantly, the model’s parameters correspond to measurable patient-specific factors, allowing us to explore how individual patient characteristics influence long-term treatment success. The parametrization of the model maps on to measurable patient characteristics, and we discuss how these parameters impact the long term behavior of the model.

Introduction: Glioblastoma tumors with a high migration to proliferation ratio (D/k ratio) are more drug resistant, suggesting 1) that estimating D/k ratios can help predict patient responses, and 2) investigating the cellular mechanisms behind D/k ratios can help understand the mechanisms of drug resistance. Here we quantify the identifiability of D/k ratios using synthetic tumor data. Materials and Methods: We used the standard reaction-diffusion model (i.e. logistic growth and diffusion) initialized with a single cell. Using this model, we grew synthetic tumors, saving the density field for calibration once the tumor had filled 10% of the domain. We fixed proliferation (k, 1/day) and used grid search followed by Levenberg-Marquardt optimization to calibrate the diffusion rate (D, mm2/day) and growth time (t, days) that produced density distributions matching the synthetic data as closely as possible. We quantified the ability of the algorithm to accurately and precisely identify D/k ratios in both the presence and absence of Gaussian noise. Results: Our algorithm finds D/k ratios with very high accuracy: 0.95 +/- 0.72% difference from correct D/k with k fixed at 0.01/day. We found that with 5% noise added the ability to accurately recover D/k ratios improved as its magnitude increased: 9.6 +/- 13.38% when D/k = 10-2 mm2 versus 2.7 +/- 3.65% when D/k = 1 mm2 Future Directions: Ongoing work will establish the minimal data requirements to accurately estimate D/k ratios within 10% of the correct values, and then apply the technique to the brain tumor image segmentation (BRATS) dataset. All BRATS patients have a single segmented pre-treatment MRI (N = 103) which we will use to estimate their D/k ratios. A subset of the BRATS patients also has RNA microarray data available (N = 91). We plan to correlate the patients’ D/k ratios with gene set scores derived from their microarray data to identify cellular mechanisms that potentially underly the D/k ratios.

Early prediction of response to therapy or lack thereof is essential for efficient treatment planning. Next-gen sequencing (NGS) has made it possible to non-invasively collect and sequence circulating tumor DNA (ctDNA) from longitudinal plasma samples. While ctDNA data continues to be collected on a wide variety of cancers and treatment types, it is still unclear what ctDNA biomarkers are most indicative of treatment success or failure. We present a pharmacokinetic model of ctDNA shedding under targeted therapy. Using this model to simulate a cohort of virtual patients, we demonstrate the predictive potential of a biomarker based on early ctDNA sampling. We compare the performance of the biomarker to that of a neural network classifier as well as existing ctDNA biomarkers. We show that our biomarker is able to match and exceed the performance of alternatives, both in terms of accuracy and earliness of prediction.

Paediatric diffuse midline gliomas are highly aggressive, incurable, childhood tumours. Their location in the brainstem limits treatment to radiotherapy, which allows an average survival time of 9-11 months. A sub-population of quiescent tumour cells are thought to be responsible for the poor outcomes of these patient. Quiescence is often viewed as a reversible resting state in which cells temporarily exit the cell cycle, the process controlling DNA replication and cell division. Radiosensitivity varies during the cell cycle, and quiescent cells exhibit a higher relative level of radio-resistance. Hypoxia (physiologically low levels of oxygen) also impacts cell cycle progression and quiescence, as well as response to radiotherapy, contributing to poor patient outcomes. We build on existing mathematical models of cell cycle progression under treatment which account for the radio-resistance of quiescent cells and their ability to re-enter the cell cycle and proliferate. We derive a system of partial differential equations (PDEs), which structures cells by the time spent in each cell cycle phase and allows transitions to and from a quiescent phase. By considering oxygen-dependent cell cycle progression, we use the model to investigate how the proportion of quiescent cells changes when we impose fluctuating oxygen dynamics and treat with radiotherapy. We extend existing studies that optimise treatment schedules using a balance of treatment outcome and toxicity/cost by incorporating a fixed radiotherapy schedule to investigate the impact of a hypothetical drug that alters the transition dynamics to and/or from quiescence. By considering different mechanisms of action for this hypothetical drug, we use our PDE model to identify candidate drugs with the potential to slow tumour progression and improve patient outcomes. This work will inform our clinical collaborators if such an improvement is possible, and what the design of a suitable drug should be.

Concurrent chemoradiation followed by brachytherapy is the standard-of-care treatment for locally advanced cervical cancer (LACC), but 30% of treated patients experience local recurrence [1], indicating a need for patient-specific, optimized therapeutic regimens to improve outcomes. We aim to predict patient-specific response to chemoradiation by applying our established MRI-based mathematical modeling framework [2]. The study cohort consisted of 10 LACC patients who underwent imaging with T2-weighted MRI, dynamic contrast-enhanced MRI, and diffusion-weighted MRI (DWI) before (V1), after two weeks (V2), and after five weeks (V3) of chemoradiation [3]. We registered all patient-specific MRI data within and between visits, and maps of the number of tumor cells (NTC) were calculated from the DWI-derived apparent diffusion coefficients. Our biology-based reaction-diffusion models characterize the spatiotemporal change in NTC as a function of cell diffusion, proliferation, and therapy-induced death. We consider two model options for cell death: (A) distinct chemotherapy (exponential decay) and radiotherapy (instantaneous decrease according to the linear-quadratic model) terms and (B) a single exponential decay term describing chemoradiation. Proliferation and chemoradiation efficacy rates were calibrated to the V1 and V2 NTC maps, and the calibrated model was run forward to predict the NTC at V3. Using Model (A), the concordance correlation coefficient (CCC) between the observed and predicted V1 to V3 change in total tumor cellularity was 0.87, and the CCC between the observed and predicted change in tumor volume was 0.90; using Model (B), the CCC values were 0.97 and 0.90, respectively. These preliminary findings show the promise of our mathematical modeling framework in predicting LACC response to chemoradiation. References: [1]. CCCMAC. Cochrane Database Syst Rev. 2010. [2]. Jarrett et al. Nat Protoc. 2021. [3]. Bowen et al. J Magn Reson Imaging. 2018.

Effective cancer therapies, while continuously improving, are often constrained by lower detection limits of disease. Tumor-immune dynamics in this limit present one of the most pressing knowledge gaps as cancer ultimate escape or elimination are often determined following an intervening period of population equilibrium sustained at low population size. Population dynamics in this small-population limit are affected by intrinsic noise in the tumor-immune interaction, as are estimates of population disease burden by extrinsic noise in acquiring such estimates through associated biomarkers. We present a modeling framework that investigates the interactions between tumor cells and the immune system in the small population regime, focusing on how these interactions influence biomarker levels. The framework combines deterministic elements, which describe tumor growth and immune responses, with stochastic components that capture the inherent variability in biomarker release. We use a system of ordinary differential equations (ODEs) to represent the tumor-immune dynamics between an adaptive immune compartment, immunogenic tumor cells, and evasive tumor cells. The immune system’s role in controlling tumor growth is reflected in the tumor-immune interaction terms. Apoptotic death via tumor-immune interactions and necrotic death via the tumor competition under a shared carrying capacity both contribute to the release of a tumor biomarker. We focus on applying our model to ct-DNA, wherein we frame ct-DNA dynamics using a stochastic differential equation (SDE). This SDE framework accounts for the variability in ct-DNA release due to the dynamic tumor-immune interactions, as well as inherent biological noise, such as DNA degradation and clearance. By coupling the ODE system of equations for tumor-immune dynamics with the SDE for ct-DNA release, we can use the model to study the fluctuations in ct-DNA levels driven by tumor-immune dynamics and exogenous sampling noise.

Melanoma cells exhibit a continuum of proliferative to invasive phenotypes. While single-cell and spatial transcriptomics have enabled biologists to quantify the distribution of phenotype amongst melanoma cells, a complete mechanistic understanding is currently lacking. A key issue is the impact of adipocyte-derived lipids, whose uptake by melanoma cells drives an invasive response that may lead to metastasis. To address this, we have developed a phenotype-structured model for melanoma cell populations that couples the phenotype dynamics to the essential aspects of intracellular lipid metabolism and the extracellular microenvironment. In this talk, I will first introduce a single-cell ODE model that illustrates how lipid uptake gives rise to phenotypic bistability in melanoma cells. I will then show how a phenotype-structured population model, whose advection term is informed by the single-cell model, exhibits a range of qualitative behaviours, including cyclic solutions and bimodal phenotypic distributions. Together, these results increase understanding of the role played by adipocyte-derived lipids and other microenvironment factors in shaping the distribution of phenotype in melanoma cell populations. We speculate that our modelling framework may also be applicable to other lipid-rich tumours (e.g. breast and ovarian cancers) that are commonly associated with increased metastasis.

Metabolite transport is essential for cellular homeostasis, energy production, and metabolic adaptation. In cancer, dysregulated transport sustains tumor growth and alters redox balance. The mitochondrial solute carrier SLC25A10 facilitates succinate, malate, and phosphate exchange, influencing central carbon metabolism. However, its transport kinetics and physiological directionality remain poorly understood. We present a mathematical model of SLC25A10 based on a ping-pong kinetic mechanism, capturing competitive dynamics between malate and succinate. Our simulations reveal that under normal conditions, malate flux dominates due to its higher binding affinity. However, in succinate dehydrogenase (SDH) dysfunction, excess succinate induces a transient efflux shift and phosphate flux reversal. If experimentally validated, this metabolic shift could serve as a biomarker for tumors with SDH mutations. Integrating our kinetic model with genome-scale metabolic networks, we highlight the role of mitochondrial transport in cancer metabolism. Specifically, in multiple myeloma, metabolic crosstalk between plasma cells and bone marrow stromal cells is key to tumor progression. Our findings demonstrate the power of mathematical modeling in uncovering transport-mediated metabolic vulnerabilities, offering potential therapeutic targets for cancer and metabolic diseases.

Blood cell formation is a complex process which is driven by hematopoietic stem cells (HSCs). HSCs give rise to progenitors and precursors which eventually produce mature blood cells, such as white blood cells, red blood cells, and platelets. Acute myeloid leukemia (AML) is an aggressive blood cancer which originates from leukemic stem cells (LSCs) and is characterized by the accumulation of aberrant immature cells, referred to as leukemic blasts. Due to the impairment of healthy blood cell formation, many AML patients suffer from life-threatening complications, such as bleeding or infection. Although treated with high-dose chemotherapy, many patients relapse and need salvage therapy. To reveal the mechanisms of disease progression and relapse, we proposed a mathematical model that accounts for competition of HSCs and LSCs in the stem cell niche and physiological feedback regulations before, during, and after chemotherapy. We fit the model to data of 7 individual patients and simulate variations of the treatment protocol. Our simulation results can recapitulate the non-monotonic recovery of HSCs observed in relapsing patients. The model suggests using the decline of HSC counts during remission as an indication for salvage therapy in patients lacking minimal residual disease markers. To bring our model closer to clinical applications, we propose a machine learning assisted mechanistic model that ensuring adherence to biological principles while learning from a larger clinical AML dataset. By embedding mechanistic constraints into machine learning, we aim to identify patient-specific predictors of relapse while preserving biological interpretability.

MD Anderson's Enderling lab recently invented a spectral spatial analysis method for estimating tumor cell diffusivity and proliferation rate from single-point-in-time biopsies of breast cancer. In combination with clinical data from the patients these parameters could help identify a new biomarker for radiotherapy. In their first study, they investigate the relationship between the power spectral density (PSD) of the three-dimensional reaction-diffusion (RD) equation with exponential growth (as model of spreading cancer cells) and the two-point correlation function of the cell distribution in the biopsy (a spatial statistic). Their results make the approach seem promising, and this work aims to validate and extend their findings. Firstly, we develop a model to generate in-silico data to validate the parameter estimation method. This is done by solving the RD equation for different growth terms (exponential and logistic), adding Gaussian noise and 'translating' its continuous results into spatial point patterns which are interpreted as cell nuclei in the 'biopsy', and then applying the method to see if the original parameters can be retrieved. This model contains several features: dimensionality can be switched between 2D and 3D, cell size can be adjusted, cuts can be added to the point pattern, and in the 3D case, biopsy thickness is variable and the plane where the slice through the 'tumor' is made can be freely chosen. And secondly, the spectral analysis method is altered by proposing a numerical solution to the PSD of the RD equation with logistic growth (valid for arbitrary dimensions). Logistic growth is assumed to be the more realistic model, however, it is harder to handle as no analytical solution is available for the equation, and hence neither for the PSD. The validation results from the in-silico data are assessed and their meaning for the application to real patient data is discussed under consideration of the different types of cell growth.

Tertiary lymphoid structures (TLSs) are organised aggregates of immune cells that form at sites of inflammation in chronic diseases, such as cancer. It is hypothesised that, in cancer, TLSs act as local hubs for the generation and regulation of a tumour-specific immune response from inside the tumour microenvironment (TME). TLSs initially form as well-mixed aggregates of T- and B-cells and mature into organised structures consisting of an inner B-cell zone surrounded by an outer T-cell zone. The presence of TLSs correlates with positive patient outcomes in several cancer types, but the mechanisms governing their formation, maturation, and role in the antitumour response remain poorly understood. Motivated by analysis of spatial transcriptomics images of TLSs in colorectal cancer, we develop an agent-based model to investigate TLS formation, maturation, and function in cancer. We model T-cells and B-cells as discrete agents which are attracted to diffusible chemokines (CXCL13 and CCL19) produced by resident stromal cells in the TME. These interactions lead to the formation of a well-mixed lymphoid aggregate that later matures into distinct T- and B-cell zones due to the segregated expression of these chemokines. Our results identify key parameters governing TLS development and suggest conditions under which TLSs are able to control tumour growth. This framework provides a qualitative basis for understanding TLS dynamics and their potential role in cancer immunotherapy.

Effective cancer therapies, while continuously improving, are often constrained by lower detection limits of disease. Tumor-immune dynamics in this limit present one of the most pressing knowledge gaps as cancer ultimate escape or elimination are often determined following an intervening period of population equilibrium sustained at low population size. Population dynamics in this small-population limit are affected by intrinsic noise in the tumor-immune interaction, as are estimates of population disease burden by extrinsic noise in acquiring such estimates through associated biomarkers. We present a modeling framework that investigates the interactions between tumor cells and the immune system in the small population regime, focusing on how these interactions influence biomarker levels. The framework combines deterministic elements, which describe tumor growth and immune responses, with stochastic components that capture the inherent variability in biomarker release. We use a system of ordinary differential equations (ODEs) to represent the tumor-immune dynamics between an adaptive immune compartment, immunogenic tumor cells, and evasive tumor cells. The immune system’s role in controlling tumor growth is reflected in the tumor-immune interaction terms. Apoptotic death via tumor-immune interactions and necrotic death via the tumor competition under a shared carrying capacity both contribute to the release of a tumor biomarker. We focus on applying our model to ct-DNA, wherein we frame ct-DNA dynamics using a stochastic differential equation (SDE). This SDE framework accounts for the variability in ct-DNA release due to the dynamic tumor-immune interactions, as well as inherent biological noise, such as DNA degradation and clearance. By coupling the ODE system of equations for tumor-immune dynamics with the SDE for ct-DNA release, we can use the model to study the fluctuations in ct-DNA levels driven by tumor-immune dynamics and exogenous sampling noise.

Melanoma cells exhibit a continuum of proliferative to invasive phenotypes. While single-cell and spatial transcriptomics have enabled biologists to quantify the distribution of phenotype amongst melanoma cells, a complete mechanistic understanding is currently lacking. A key issue is the impact of adipocyte-derived lipids, whose uptake by melanoma cells drives an invasive response that may lead to metastasis. To address this, we have developed a phenotype-structured model for melanoma cell populations that couples the phenotype dynamics to the essential aspects of intracellular lipid metabolism and the extracellular microenvironment. In this talk, I will first introduce a single-cell ODE model that illustrates how lipid uptake gives rise to phenotypic bistability in melanoma cells. I will then show how a phenotype-structured population model, whose advection term is informed by the single-cell model, exhibits a range of qualitative behaviours, including cyclic solutions and bimodal phenotypic distributions. Together, these results increase understanding of the role played by adipocyte-derived lipids and other microenvironment factors in shaping the distribution of phenotype in melanoma cell populations. We speculate that our modelling framework may also be applicable to other lipid-rich tumours (e.g. breast and ovarian cancers) that are commonly associated with increased metastasis.

Metabolite transport is essential for cellular homeostasis, energy production, and metabolic adaptation. In cancer, dysregulated transport sustains tumor growth and alters redox balance. The mitochondrial solute carrier SLC25A10 facilitates succinate, malate, and phosphate exchange, influencing central carbon metabolism. However, its transport kinetics and physiological directionality remain poorly understood. We present a mathematical model of SLC25A10 based on a ping-pong kinetic mechanism, capturing competitive dynamics between malate and succinate. Our simulations reveal that under normal conditions, malate flux dominates due to its higher binding affinity. However, in succinate dehydrogenase (SDH) dysfunction, excess succinate induces a transient efflux shift and phosphate flux reversal. If experimentally validated, this metabolic shift could serve as a biomarker for tumors with SDH mutations. Integrating our kinetic model with genome-scale metabolic networks, we highlight the role of mitochondrial transport in cancer metabolism. Specifically, in multiple myeloma, metabolic crosstalk between plasma cells and bone marrow stromal cells is key to tumor progression. Our findings demonstrate the power of mathematical modeling in uncovering transport-mediated metabolic vulnerabilities, offering potential therapeutic targets for cancer and metabolic diseases.

Blood cell formation is a complex process which is driven by hematopoietic stem cells (HSCs). HSCs give rise to progenitors and precursors which eventually produce mature blood cells, such as white blood cells, red blood cells, and platelets. Acute myeloid leukemia (AML) is an aggressive blood cancer which originates from leukemic stem cells (LSCs) and is characterized by the accumulation of aberrant immature cells, referred to as leukemic blasts. Due to the impairment of healthy blood cell formation, many AML patients suffer from life-threatening complications, such as bleeding or infection. Although treated with high-dose chemotherapy, many patients relapse and need salvage therapy. To reveal the mechanisms of disease progression and relapse, we proposed a mathematical model that accounts for competition of HSCs and LSCs in the stem cell niche and physiological feedback regulations before, during, and after chemotherapy. We fit the model to data of 7 individual patients and simulate variations of the treatment protocol. Our simulation results can recapitulate the non-monotonic recovery of HSCs observed in relapsing patients. The model suggests using the decline of HSC counts during remission as an indication for salvage therapy in patients lacking minimal residual disease markers. To bring our model closer to clinical applications, we propose a machine learning assisted mechanistic model that ensuring adherence to biological principles while learning from a larger clinical AML dataset. By embedding mechanistic constraints into machine learning, we aim to identify patient-specific predictors of relapse while preserving biological interpretability.

MD Anderson's Enderling lab recently invented a spectral spatial analysis method for estimating tumor cell diffusivity and proliferation rate from single-point-in-time biopsies of breast cancer. In combination with clinical data from the patients these parameters could help identify a new biomarker for radiotherapy. In their first study, they investigate the relationship between the power spectral density (PSD) of the three-dimensional reaction-diffusion (RD) equation with exponential growth (as model of spreading cancer cells) and the two-point correlation function of the cell distribution in the biopsy (a spatial statistic). Their results make the approach seem promising, and this work aims to validate and extend their findings. Firstly, we develop a model to generate in-silico data to validate the parameter estimation method. This is done by solving the RD equation for different growth terms (exponential and logistic), adding Gaussian noise and 'translating' its continuous results into spatial point patterns which are interpreted as cell nuclei in the 'biopsy', and then applying the method to see if the original parameters can be retrieved. This model contains several features: dimensionality can be switched between 2D and 3D, cell size can be adjusted, cuts can be added to the point pattern, and in the 3D case, biopsy thickness is variable and the plane where the slice through the 'tumor' is made can be freely chosen. And secondly, the spectral analysis method is altered by proposing a numerical solution to the PSD of the RD equation with logistic growth (valid for arbitrary dimensions). Logistic growth is assumed to be the more realistic model, however, it is harder to handle as no analytical solution is available for the equation, and hence neither for the PSD. The validation results from the in-silico data are assessed and their meaning for the application to real patient data is discussed under consideration of the different types of cell growth.

Tertiary lymphoid structures (TLSs) are organised aggregates of immune cells that form at sites of inflammation in chronic diseases, such as cancer. It is hypothesised that, in cancer, TLSs act as local hubs for the generation and regulation of a tumour-specific immune response from inside the tumour microenvironment (TME). TLSs initially form as well-mixed aggregates of T- and B-cells and mature into organised structures consisting of an inner B-cell zone surrounded by an outer T-cell zone. The presence of TLSs correlates with positive patient outcomes in several cancer types, but the mechanisms governing their formation, maturation, and role in the antitumour response remain poorly understood. Motivated by analysis of spatial transcriptomics images of TLSs in colorectal cancer, we develop an agent-based model to investigate TLS formation, maturation, and function in cancer. We model T-cells and B-cells as discrete agents which are attracted to diffusible chemokines (CXCL13 and CCL19) produced by resident stromal cells in the TME. These interactions lead to the formation of a well-mixed lymphoid aggregate that later matures into distinct T- and B-cell zones due to the segregated expression of these chemokines. Our results identify key parameters governing TLS development and suggest conditions under which TLSs are able to control tumour growth. This framework provides a qualitative basis for understanding TLS dynamics and their potential role in cancer immunotherapy.

The complexity of mathematical models describing respiratory mechanics has grown in recent years, however, parameter identifiability of such models has only been studied in the last decade in the context of observable data. This study investigates parameter identifiability of a nonlinear respiratory mechanics model tuned to the physiology of an extremely preterm infant, using global Morris screening, local deterministic sensitivity analysis, and singular value decomposition-based subset selection. The model predicts airflow and dynamic pulmonary volumes and pressures under varying levels of continuous positive airway pressure, and a range of parameters characterizing both surfactant-treated and surfactant-deficient lung. The model was adapted to data from a spontaneously breathing 1 kg infant using gradient-based optimization to estimate the parameter subset characterizing the patient's state of health.

As consumer-grade wearable technology has become more prevalent in recent years, large-scale collections of data have been made available for researchers. We analyze significant amounts of wearable data to determine the circadian features that differ across groups and time frames. Using wearable activity, steps, and heart rate data, we adapt mathematical models to accurately estimate circadian phase across populations. This has a number of applications, including chronotherapeutic drug delivery, reducing fatigue, and shift work scheduling. We demonstrate applications to estimating circadian phase (dim light melatonin onset, or DLMO) in a home-based cohort of later-life adults, showing that activity-based models perform similarly or better than light-based models in DLMO estimation. We further use these models to provide wearable-based lighting interventions for reducing cancer-related fatigue. In particular, we test whether these lighting interventions, delivered via a mobile app, reduce cancer-related fatigue in a randomized controlled trial with 138 breast cancer, prostate cancer, and hematopoietic stem cell transplant patients. These interventions, based on real-time assessment of circadian rhythm through wearable devices, improve certain measures of fatigue (e.g., daily measurements of fatigue) in cancer patients. Further studies are needed to tune these models and assess the effect of lighting interventions in broader and more diverse cancer care settings.

Fetal weight estimation via ultrasound is performed by measuring biometric parameters such as the biparietal diameter (BPD), head circumference (HC), abdominal circumference (AC), and femur length (FL), which are then used in various mathematical formulas to calculate the estimated weight. But do all parameters matter equally? To assess their contribution on fetal weight estimation, we analyzed 29 published formulas across 26 diverse global datasets. Results show that AC is consistently the parameter of greatest importance, while head measurements (BPD, HC) often add little value, particularly in the later stages of pregnancy. Additionally, nearly half of the formulas include redundant parameters, and two-thirds exhibit a crossover in parameter importance—some transition from low to high significance, while others decline from high to low—over the course of gestation. These findings highlight opportunities to simplify fetal weight estimation for clinicians, prioritizing AC reliability and trimming unnecessary inputs. Our work bridges mathematics and prenatal care, offering clearer guidelines to improve ultrasound-based predictions and support healthier pregnancy outcomes.