Atherosclerosis is a chronic inflammatory disease characterised by lipid accumulation within arterial walls and driven by macrophage interactions with extracellular material, particularly lipid. In this work, we use mathematical modelling to investigate the dynamics of the macrophages in early atherosclerosis. We develop a discrete, lipid-structured mathematical model that accounts for cell proliferation and crowding, and also extracellular material uptake and offloading. With this model we are able to describe the dynamics of key biophysical quantities in the plaque, in particular the total number of macrophages and the total amount of intracellular material contained within the macrophages. Moreover, we rigorously derive a continuum approximation of the discrete model using the method of discrete multiple scales and asymptotic analysis techniques. In this way, we systematically derive a partial differential equation that accurately describes the distribution of macrophage content at leading order. We take advantage of the continuum form to analyse the mathematical system and understand its biological implications, such as the effects of proliferation and crowding on plaque composition. We also investigate the important spatio-temporal regions that appear degenerate but can be understood via boundary layer analysis.
Atherosclerosis is a chronic inflammatory disease characterised by lipid accumulation within arterial walls and driven by macrophage interactions with extracellular material, particularly lipid. In this work, we use mathematical modelling to investigate the dynamics of the macrophages in early atherosclerosis. We develop a discrete, lipid-structured mathematical model that accounts for cell proliferation and crowding, and also extracellular material uptake and offloading. With this model we are able to describe the dynamics of key biophysical quantities in the plaque, in particular the total number of macrophages and the total amount of intracellular material contained within the macrophages. Moreover, we rigorously derive a continuum approximation of the discrete model using the method of discrete multiple scales and asymptotic analysis techniques. In this way, we systematically derive a partial differential equation that accurately describes the distribution of macrophage content at leading order. We take advantage of the continuum form to analyse the mathematical system and understand its biological implications, such as the effects of proliferation and crowding on plaque composition. We also investigate the important spatio-temporal regions that appear degenerate but can be understood via boundary layer analysis.
Understanding tumour dynamics under hypoxic conditions is critical for optimising cancer therapies, particularly with chemotherapeutic agents like Paclitaxel. This study presents a refined mathematical model of tumour growth that incorporates Paclitaxel effects and hypoxia-driven resistance using a system of nonlinear ordinary differential equations (ODEs). We employ the Metropolis-Hastings Markov Chain Monte Carlo (MH MCMC) algorithm for Bayesian inversion and parameter estimation, providing a probabilistic framework to capture uncertainties. Sensitivity analysis is conducted using the multiple shooting method, which enhances the stability and accuracy of local sensitivity estimates across time intervals. The simulation results demonstrate that cell viability is reduced under moderate hypoxia when treated with Paclitaxel, which is consistent with experimental data from HCC1806 breast cancer cell lines. This agreement between model predictions and experimental outcomes supports the model’s validity in capturing key biological mechanisms. Future work will extend the model using Physics-Informed Neural Networks (PINNs) to improve computational efficiency and explore advanced inverse problem-solving techniques for robust cancer treatment optimisation.
Spatial gene-expression data can be clustered to segment a tissue into distinct spatial domains representing tissue structure. Though clustering algorithms are limited to a single fixed scale (by choice of a resolution hyperparameter k), we develop new methods from topological data analysis to analyze patterns in clusters across multiple scales. Zero-dimensional persistent homology analyzes the connectivity of data by tracking changes in homology groups across a filtered simplicial complex. We build a new filtration scheme to analyze similarity between clusters generated from multiple choices of scale parameter k, where persistent components represent clusters which exist across scales. We apply these results to select optimal scale parameters for spatial gene-expression clustering. These results have potential clinical application in tumor identification, where the size and scale of cancerous domains within healthy tissue is not known a priori.
Abstract A biochemical system includes a network of chemical reactions often exhibiting complex behaviors such as oscillations, spatial patterns, and multistability. The parameter values of these models are often unknown or difficult to measure, and even some details of the reaction networks may be uncertain. Since these models tend to be large and complex, it is useful to create a simplified version of these models. However, traditional model-reduction methods depend on knowledge of parameter values which make them difficult to apply. Qualitative stability analysis methods provide an alternative approach without necessarily requiring parameter values. When reducing models with non-trivial dynamics arising from an instability, one must ensure that the conditions for instability are preserved, which depend mainly on the presence of circuits, and their signs. Roussel and Soares presented dynamics-preserving reductions based on Ivanova's qualitative conditions for instabilities (J. Math. Biol. 89, 42). The main objective of this research is to implement a similar framework based on the concepts outlined in that paper. However, instead of using Ivanova's conditions for instability, we will apply the Thomas qualitative stability analysis method to preserve the structures in the interaction graph that generate instability. An Oregonator-class model for oscillations in the photosensitive Belousov-Zhabotinsky (BZ) reaction due to Amemiya and coworkers is used in an initial exploration of possible reduction rules in interaction graphs. Given that the interaction graph discards information about the kinetics of a reaction, some attention will have to be given to the potential loss of important nonlinear terms while implementing the new method.
Repeated patterns such as bristles and hair follicles play an important role in epithelia, which sense the environment. Optimal organization of patterns contributes to normal tissue function and gives organisms a spatial and temporal mapped-out input of their environmental stimuli. Although many local (e.g., cell- cell) signaling mechanisms are understood, some gaps still exist in our understanding of long-distance signaling via cell protrusions such as filopodia and cytoneme. The sensory bristles of the fruit fly Drosophila Melanogaster are a genetically tractable system for studying the formation of repeating patterns and invariably long-range cell signaling via cell protrusions. One critical feature of the sensory bristle spot pattern is the presence of long-range lateral inhibition, a mechanism that relies on forming actin-based cell protrusions – filopodia. We develop a mathematical model to describe filopodia dynamics and their role in determining cell fate during patterning.
The endoplasmic reticulum (ER) is a highly dynamic organelle that undergoes contin- uous remodeling between tubular and sheet-like structures, driven by the Rtn and Reep protein families. Understanding the physical principles underlying these transitions is cru- cial for elucidating the ER’s role in cellular homeostasis and disease. In this study, we em- ploy mesoscale simulations to investigate the mechanisms by which curvature-promoting proteins regulate ER morphology. Specifically, we explore the influence of protein in- trinsic curvature, protein concentration, and protein sti!ening on tubulation dynamics. Our results indicate that increasing the intrinsic curvature of proteins lowers the pro- tein coverage threshold required for tubulation, while enhanced membrane sti!ness facil- itates curvature propagation at lower protein coverage. A phase diagram is constructed to map the conditions necessary for membrane remodeling, identifying critical curvature and protein coverage thresholds that drive ER transformation. These findings establish a quantitative framework for ER shape regulation, shedding light on the interplay between protein-membrane interactions and mechanical properties in ER morphogenesis. By inte- grating computational predictions, this study advances our understanding of ER structural dynamics and its implications for cellular function.
Antibiotics, which can be defined as substances that work against bacteria, are one of the most useful agents used in healthcare. As a result, they serve to treat and prevent many bacterial infections. However, due to the emergence of antibiotic resistance, where bacteria develop a mechanism to defend themselves against antibiotics, managing infections has become increasingly challenging. Antibiotic resistance in bacteria arises through genetic mutations or horizontal gene transfer. Spatial heterogeneity in antibiotic concentration has a potential to affect this bacterial evolution. For example, compared to a well mixed population, in a highly structured population, increased phenotypic and genotypic diversity, as well as slower adaptation, is expected. Here, we are studying the bacterial evolution under the stochastic processes of division, which is influenced by the availability of food sources in the culture, as well as by mutations and migration. As division reaction is time dependent, this chemical system is non-homogeneous and non-stationary. In this scenario, continuous time Markov processes can not be applied as chemical reactions are non-homogeneous and non-stationary. In this study, an expression was formulated to determine the time until the next reaction occurs, given the current state of the system, by considering the combined effects of division, migration, and mutation.
Replication stress refers to the impeding of DNA copying and the slowing or arresting of replication forks during DNA synthesis. It arises due to a number of exogenous and endogenous agents such as reactive oxygen species (ROS), radiation-induced DNA lesions, and noncanonical folded DNA species such as G-quadruplexes. Replication stress can give rise to chromosomal missegregation in anaphase, DNA breakage, or faulty rearrangements. In this work, we take a stochastic process approach to modeling replication stress, using a coupled system of point processes to capture replication fork distribution and characterization of replication origin licensing, and Poisson process modeling of origin activation. Recent work in yeast has demonstrated the appropriateness of the Poisson model for capturing the stochastic multiple activation process under replication stress. Finally we focus particularly on the role of G-quadruplexes (G4), which are guanine (G)-rich regions of DNA that form noncanonical quadruple-stranded structures that are implicated in replication stress. We discuss how the different topologies of G4 potentially influence the origin activation process modeled in this work.
Reaction-diffusion systems are fundamental in modeling the complex spatiotemporal dynamics in biological, chemical, and ecological phenomena. In this study, we investigate a bistable reaction-diffusion system motivated by the experimental observations on Rho-GEF-Myosin signaling network that controls cell contraction dynamics. Through a combination of numerical bifurcation analysis and simulations, we explore how diffusion alters the intrinsic dynamics of distinct temporal regimes exhibited by the underlying reaction kinetics. Our results demonstrate that: (i) diffusion can destabilize a uniform stable steady state, leading to classical Turing patterns; (ii) in oscillatory regimes, diffusion drives the system away from temporal periodicity into spatially heterogeneous oscillations, indicating far-from-equilibrium behavior; and (iii) in bistable regions, diffusion induces pattern formation, wave propagation, and oscillatory pulses. Floquet theory is used to quantify the diffusion-driven destabilization of a homogeneously stable limit cycle, identifying critical diffusion coefficients for diffusion-driven instability. These findings offer theoretical insights into diffusion-induced transitions and can contribute to the broader understanding of pattern formation and dynamic regulation in developmental and cellular biology.
Reaction-diffusion systems are fundamental in modeling the complex spatiotemporal dynamics in biological, chemical, and ecological phenomena. In this study, we investigate a bistable reaction-diffusion system motivated by the experimental observations on Rho-GEF-Myosin signaling network that controls cell contraction dynamics. Through a combination of numerical bifurcation analysis and simulations, we explore how diffusion alters the intrinsic dynamics of distinct temporal regimes exhibited by the underlying reaction kinetics. Our results demonstrate that: (i) diffusion can destabilize a uniform stable steady state, leading to classical Turing patterns; (ii) in oscillatory regimes, diffusion drives the system away from temporal periodicity into spatially heterogeneous oscillations, indicating far-from-equilibrium behavior; and (iii) in bistable regions, diffusion induces pattern formation, wave propagation, and oscillatory pulses. Floquet theory is used to quantify the diffusion-driven destabilization of a homogeneously stable limit cycle, identifying critical diffusion coefficients for diffusion-driven instability. These findings offer theoretical insights into diffusion-induced transitions and can contribute to the broader understanding of pattern formation and dynamic regulation in developmental and cellular biology.
Understanding tumour dynamics under hypoxic conditions is critical for optimising cancer therapies, particularly with chemotherapeutic agents like Paclitaxel. This study presents a refined mathematical model of tumour growth that incorporates Paclitaxel effects and hypoxia-driven resistance using a system of nonlinear ordinary differential equations (ODEs). We employ the Metropolis-Hastings Markov Chain Monte Carlo (MH MCMC) algorithm for Bayesian inversion and parameter estimation, providing a probabilistic framework to capture uncertainties. Sensitivity analysis is conducted using the multiple shooting method, which enhances the stability and accuracy of local sensitivity estimates across time intervals. The simulation results demonstrate that cell viability is reduced under moderate hypoxia when treated with Paclitaxel, which is consistent with experimental data from HCC1806 breast cancer cell lines. This agreement between model predictions and experimental outcomes supports the model’s validity in capturing key biological mechanisms. Future work will extend the model using Physics-Informed Neural Networks (PINNs) to improve computational efficiency and explore advanced inverse problem-solving techniques for robust cancer treatment optimisation.
Osteocytes are specialized bone cells responsible for sensing mechanical cues and directing bone remodeling. These cells reside in small cavities called lacunae within the bone matrix and extend dendritic processes through narrow channels known as canaliculi to connect with other osteocytes. Surrounding the cell and lining the lacuno-canalicular system is an interstitial fluid and cellular coating called the pericellular matrix (PCM). Previous studies have shown that the stress and strain required to elicit a significant response in osteocytes are approximately ten times greater than those typically experienced during normal physical activity. However, the mechanism by which macroscale mechanical signals are amplified to such levels within the osteocyte network is not yet fully understood. We develop coarse-grained models to investigate fluid-osteocyte interactions. These models incorporate key components of the osteocyte and its microenvironment, including the cell body (comprising the membrane/cortex, cytoskeleton, and cytosol), cellular processes, canaliculi, lacuna, interstitial fluid, and the surrounding bone matrix. The cell membrane is represented by a cross-linked viscoelastic fiber network arranged in triangular patterns. The cellular processes are modeled by discretized curves using damped elastic springs, while the cytoskeleton is constructed from tetrahedral elements, with each edge represented by a linearly viscoelastic fiber. Both the cytosol and interstitial fluid are treated as viscous, incompressible fluids. The surrounding bone is modeled as a rigid material. Both extracellular and intracellular flows are governed by the Navier-Stokes equations and numerically solved by the lattice Boltzmann method. The fluid-osteocyte interactions are simulated using the immersed boundary method. A key finding from our simulations is that stress and strain are highly concentrated at the junctions where the cellular processes connect to the main body of the cell.
Spatial gene-expression data can be clustered to segment a tissue into distinct spatial domains representing tissue structure. Though clustering algorithms are limited to a single fixed scale (by choice of a resolution hyperparameter k), we develop new methods from topological data analysis to analyze patterns in clusters across multiple scales. Zero-dimensional persistent homology analyzes the connectivity of data by tracking changes in homology groups across a filtered simplicial complex. We build a new filtration scheme to analyze similarity between clusters generated from multiple choices of scale parameter k, where persistent components represent clusters which exist across scales. We apply these results to select optimal scale parameters for spatial gene-expression clustering. These results have potential clinical application in tumor identification, where the size and scale of cancerous domains within healthy tissue is not known a priori.
Abstract A biochemical system includes a network of chemical reactions often exhibiting complex behaviors such as oscillations, spatial patterns, and multistability. The parameter values of these models are often unknown or difficult to measure, and even some details of the reaction networks may be uncertain. Since these models tend to be large and complex, it is useful to create a simplified version of these models. However, traditional model-reduction methods depend on knowledge of parameter values which make them difficult to apply. Qualitative stability analysis methods provide an alternative approach without necessarily requiring parameter values. When reducing models with non-trivial dynamics arising from an instability, one must ensure that the conditions for instability are preserved, which depend mainly on the presence of circuits, and their signs. Roussel and Soares presented dynamics-preserving reductions based on Ivanova's qualitative conditions for instabilities (J. Math. Biol. 89, 42). The main objective of this research is to implement a similar framework based on the concepts outlined in that paper. However, instead of using Ivanova's conditions for instability, we will apply the Thomas qualitative stability analysis method to preserve the structures in the interaction graph that generate instability. An Oregonator-class model for oscillations in the photosensitive Belousov-Zhabotinsky (BZ) reaction due to Amemiya and coworkers is used in an initial exploration of possible reduction rules in interaction graphs. Given that the interaction graph discards information about the kinetics of a reaction, some attention will have to be given to the potential loss of important nonlinear terms while implementing the new method.
Repeated patterns such as bristles and hair follicles play an important role in epithelia, which sense the environment. Optimal organization of patterns contributes to normal tissue function and gives organisms a spatial and temporal mapped-out input of their environmental stimuli. Although many local (e.g., cell- cell) signaling mechanisms are understood, some gaps still exist in our understanding of long-distance signaling via cell protrusions such as filopodia and cytoneme. The sensory bristles of the fruit fly Drosophila Melanogaster are a genetically tractable system for studying the formation of repeating patterns and invariably long-range cell signaling via cell protrusions. One critical feature of the sensory bristle spot pattern is the presence of long-range lateral inhibition, a mechanism that relies on forming actin-based cell protrusions – filopodia. We develop a mathematical model to describe filopodia dynamics and their role in determining cell fate during patterning.
The endoplasmic reticulum (ER) is a highly dynamic organelle that undergoes contin- uous remodeling between tubular and sheet-like structures, driven by the Rtn and Reep protein families. Understanding the physical principles underlying these transitions is cru- cial for elucidating the ER’s role in cellular homeostasis and disease. In this study, we em- ploy mesoscale simulations to investigate the mechanisms by which curvature-promoting proteins regulate ER morphology. Specifically, we explore the influence of protein in- trinsic curvature, protein concentration, and protein sti!ening on tubulation dynamics. Our results indicate that increasing the intrinsic curvature of proteins lowers the pro- tein coverage threshold required for tubulation, while enhanced membrane sti!ness facil- itates curvature propagation at lower protein coverage. A phase diagram is constructed to map the conditions necessary for membrane remodeling, identifying critical curvature and protein coverage thresholds that drive ER transformation. These findings establish a quantitative framework for ER shape regulation, shedding light on the interplay between protein-membrane interactions and mechanical properties in ER morphogenesis. By inte- grating computational predictions, this study advances our understanding of ER structural dynamics and its implications for cellular function.
Antibiotics, which can be defined as substances that work against bacteria, are one of the most useful agents used in healthcare. As a result, they serve to treat and prevent many bacterial infections. However, due to the emergence of antibiotic resistance, where bacteria develop a mechanism to defend themselves against antibiotics, managing infections has become increasingly challenging. Antibiotic resistance in bacteria arises through genetic mutations or horizontal gene transfer. Spatial heterogeneity in antibiotic concentration has a potential to affect this bacterial evolution. For example, compared to a well mixed population, in a highly structured population, increased phenotypic and genotypic diversity, as well as slower adaptation, is expected. Here, we are studying the bacterial evolution under the stochastic processes of division, which is influenced by the availability of food sources in the culture, as well as by mutations and migration. As division reaction is time dependent, this chemical system is non-homogeneous and non-stationary. In this scenario, continuous time Markov processes can not be applied as chemical reactions are non-homogeneous and non-stationary. In this study, an expression was formulated to determine the time until the next reaction occurs, given the current state of the system, by considering the combined effects of division, migration, and mutation.
Replication stress refers to the impeding of DNA copying and the slowing or arresting of replication forks during DNA synthesis. It arises due to a number of exogenous and endogenous agents such as reactive oxygen species (ROS), radiation-induced DNA lesions, and noncanonical folded DNA species such as G-quadruplexes. Replication stress can give rise to chromosomal missegregation in anaphase, DNA breakage, or faulty rearrangements. In this work, we take a stochastic process approach to modeling replication stress, using a coupled system of point processes to capture replication fork distribution and characterization of replication origin licensing, and Poisson process modeling of origin activation. Recent work in yeast has demonstrated the appropriateness of the Poisson model for capturing the stochastic multiple activation process under replication stress. Finally we focus particularly on the role of G-quadruplexes (G4), which are guanine (G)-rich regions of DNA that form noncanonical quadruple-stranded structures that are implicated in replication stress. We discuss how the different topologies of G4 potentially influence the origin activation process modeled in this work.
Reaction-diffusion systems are fundamental in modeling the complex spatiotemporal dynamics in biological, chemical, and ecological phenomena. In this study, we investigate a bistable reaction-diffusion system motivated by the experimental observations on Rho-GEF-Myosin signaling network that controls cell contraction dynamics. Through a combination of numerical bifurcation analysis and simulations, we explore how diffusion alters the intrinsic dynamics of distinct temporal regimes exhibited by the underlying reaction kinetics. Our results demonstrate that: (i) diffusion can destabilize a uniform stable steady state, leading to classical Turing patterns; (ii) in oscillatory regimes, diffusion drives the system away from temporal periodicity into spatially heterogeneous oscillations, indicating far-from-equilibrium behavior; and (iii) in bistable regions, diffusion induces pattern formation, wave propagation, and oscillatory pulses. Floquet theory is used to quantify the diffusion-driven destabilization of a homogeneously stable limit cycle, identifying critical diffusion coefficients for diffusion-driven instability. These findings offer theoretical insights into diffusion-induced transitions and can contribute to the broader understanding of pattern formation and dynamic regulation in developmental and cellular biology.
We study a fluctuation test where cell colonies grow from a single cell to a specified population size before being treated. During growth, cells may acquire resistance to treatment and pass it to offspring, with a small probability. Unlike the classical Luria–Delbrück test, we allow the resistant state to revert to a drug-sensitive state, motivated by recent research on drug tolerance in cancer and microbes. This modification does not change the central part of the Luria–Delbrück distribution, where the Landau probability density function approximation still applies. However, the right tail of the distribution deviates from the power law of the Landau distribution, with the correction factor equal to the Landau cumulative distribution function. We use singular perturbation theory and asymptotic matching to derive uniformly valid approximations and describe tail corrections for populations with different initial cell states.
As temperatures rise, photosymbiotic marine species are presented with unique challenges. Exaptasia diaphana is a highly adaptable model organism for studying these challenges. Current methods for acquiring experimental data are expensive and require sacrificing the organism. This work focuses on the development of computational tools that will reduce cost and automate the acquisition of data. A pipeline consisting of a convolutional neural network and a transformer-based model is used to accomplish this task. Given input images of aiptasia colonies, this pipeline automatically produces accurate segmentations of aiptasia that can be used for experimental data acquisition (e.g. obtaining counts, measuring oral disk size and color information, etc.).
In an information-processing investment game, such as the growth of a population of organisms in a changing environment, Kelly betting maximizes the expected log rate of growth. In this work, we show that Kelly bets are closely related to optimal single-letter codes (i.e., they can achieve the rate-distortion bound with equality). Thus, natural information processing systems with limited computational resources can achieve information-theoretically optimal performance. We show that the rate-distortion tradeoff for an investment game has a simple linear bound, and that the bound is achievable at the point where the corresponding single-letter code is optimal. This interpretation has two interesting consequences. First, we show that increasing the organism's portfolio of potential strategies can lead to optimal performance over a continuous range of channels, even if the strategy portfolio is fixed. Second, we show that increasing an organism's number of phenotypes (i.e., its number of possible behaviours in response to the environment) can lead to higher growth rate, and we give conditions under which this occurs. Examples illustrating the results in simplified biological scenarios are presented.
Bivalve Transmissible Neoplasia (BTN) is an increasingly prevalent cancer spreading among bivalve species worldwide. BTN dynamics introduce complexities not common in many other infectious diseases due to its marine environment. Therefore, little is known about its transmission dynamics and population effects. This study develops a Bayesian framework to model BTN spread within Mya arenaria populations to address key gaps in our understanding, with a focus on statistical methodology for parameter inference and model development. We use a Bayesian compartmental modeling approach to infer and refine model parameters and leverage controlled laboratory and survey data. Laboratory data provide information on cancer cell emission, disease progression and environmental factor effects on the disease. Survey data from the field includes samples from East Coast sites that are used for fitting the model to disease progression over time in its natural environment. To capture the intricacies of BTN, our model framework builds on the traditional Susceptible-Exposed-Infectious (SEI) epidemiological model by incorporating cancer particle survivability components and the environmental effects of temperature. The Bayesian approach in our model development, implemented in STAN, provides parameter inferences and quantifies uncertainty in the results amidst limited or noisy ecological data. Future work will validate the model by comparing its predictions with 2025 survey data and conducting sensitivity analyses to identify key parameters. This statistical framework not only advances our understanding of BTN, but also demonstrates the applicability of Bayesian modeling in developing complex ecological and epidemiological compartmental models.
Wetlands are characterised by the interaction of soil with seasonal or permanent water bodies and serve several crucial ecological functions including flood prevention and water filtration. Moreover, serving as a boundary between land and aquatic environments, they provide diverse ecosystems for a variety of plants, animals and microbes. However, the ability of wetlands to sequester and store carbon also secures their place as the largest natural source of methane emissions. These emissions are strongly dependent on both the relative depth of the water table to the soil and the soil temperature, with submerged warm soil providing favourable anaerobic conditions for methanogenesis by microbes and detrimental to methane consumption through oxidation. Hence, wetland methane emissions are strongly susceptible to climate change, particularly changes in rainfall and temperature. We present a mathematical model to describe the stochastic movement of the water table coupled to a simple set of ODEs describing methane production, oxidation, and emission, parameterised by this water table depth. We employ this model to predict how the inherent variations in water table depth due to the soil profile leads to changing emission profiles across individual wetlands. Moreover, by exploring the sensitivity of these emissions to rainfall and temperature changes, we demonstrate during wetland conservation efforts the need to consider how climate change will influence emissions.
In this paper, a mathematical model is proposed to explain the interaction between Eastern Hemlock Trees and the invasive species Hemlock Woolly Adelgid. A system of reaction diffusion equations is used for this modeling exercise. There are at most three (3) steady states for the system of which the coexistent state is the only stable steady states for some parameter values. The model dynamics show that the solutions exhibit traveling wave solutions. In addition, a sensitivity analysis is conducted to determine the impact of model parameters. Sensitivity analysis suggests that the mortality rate of Eastern Hemlock Trees and the predation intensity of Hemlock Woolly Adelgid drive the dynamics of the interaction. Since Eastern Hemlock trees are foundation tress that provide shelter for several species they need to be protected for as long as possible. Based on the model dynamics and sensitivity analysis, it is postulated that selective and strategic removal of these trees will help curtail their destruction.
In the present work, using the theory of dynamical system, we discuss the dynamics of a cannibalistic predator prey model in the presence of linear harvesting schemes in prey (pest) population and provision of additional food to predators (natural enemies). A detailed mathematical analysis and numerical evaluations have been presented to discuss the pest free state, coexistence of species, stability, occurrence of different bifurcations (saddle-node, transcritical, Hopf, Bogdanov-Takens) and the impact of additional food and harvesting schemes on the dynamics of the system. Interestingly, we also observe that the pest population density decreases immediately even when small amount of harvesting is implemented. Also the eradication of pest population (stable pest free state) could be achieved via variation in the additional food and implemented harvesting schemes. The individual effects of harvesting parameters on the pest density suggest that the linear harvesting scheme is more effective to control the pest population rather than constant and nonlinear harvesting schemes.
Our work focuses on modelling how predator-prey dynamics are affected by the frequency and impact of forest disturbances such as forestry and wildfires and the subsequent regeneration of the forest. Our focus is on the Goshawk- Pine Squirrel predator-prey system. The coastal subspecies of the Goshawk is listed as endangered and their preferred habitat, mature, dense coniferous forests, is commonly disturbed. Climate change is increasing the frequency and severity of wildfires as well as periods of drought and heavy rainfall, leading to new challenges for sustaining Goshawk habitat. We present a modelling approach examining the dynamics of Goshawks and Squirrels within a disturbed landscape, and determine how different types of disturbance affect persistence of both species. An additional complicating factor that we take into account is the territoriality of Squirrels and Goshawks. Ultimately, we hope that the final model will inform conservation efforts for the Goshawk and, more generally, for predator-prey systems that rely on mature forest.
This research presents two key advancements in bat conservation through mathematical modelling and ecological data integration. First, we developed a mathematical model to predict the movement of Greater Horseshoe bats, allowing us to estimate their location based on known roost positions. Second, we introduced a novel algorithm that combines bat movement data with sound recordings from static microphone detectors, enabling the identification of potential bat roost locations. This method significantly reduces the area to be searched, improving conservation efficiency and supporting sustainable building practices. Our approach not only enhances bat protection but also offers valuable tools for ecologists to optimise monitoring efforts in diverse environments.
We study a fluctuation test where cell colonies grow from a single cell to a specified population size before being treated. During growth, cells may acquire resistance to treatment and pass it to offspring, with a small probability. Unlike the classical Luria–Delbrück test, we allow the resistant state to revert to a drug-sensitive state, motivated by recent research on drug tolerance in cancer and microbes. This modification does not change the central part of the Luria–Delbrück distribution, where the Landau probability density function approximation still applies. However, the right tail of the distribution deviates from the power law of the Landau distribution, with the correction factor equal to the Landau cumulative distribution function. We use singular perturbation theory and asymptotic matching to derive uniformly valid approximations and describe tail corrections for populations with different initial cell states.
The driving mechanisms behind animal movement behaviour are expected to shift over time, for example, as the seasons progress the locations of food resources change. Often, the change points between behaviours are estimated using telemetry data. However, identifying these change points in time and corresponding drivers is challenging using telemetry data alone. Currently, several existing statistical methods are used to address this problem but are subject to strong assumptions and limitations on any inferences. We test some of the popular methods on pseudo-data with known change points and known error to compare the accuracy of the statistical models against each other. We find that existing methods perform poorly under a number of common scenarios and then success rate is highly variable. Assessing the tools in the movement behaviour toolbox, as well as expanding them, is important for the development of the field and can identify the best model for use under certain conditions.
Flowers and pollinators have co-evolved intricate strategies to find each other, which may be studied through mathematical modeling. In this talk, I present stochastic models that describe the dynamics of insect attraction to flowers. A key case study is the world’s largest flower, Rafflesia, which relies on volatile cues to lure carrion flies. Using stochastic methods, I examine how environmental variability, flower geometry and pollinator behavior influence attraction dynamics. This framework also extends to other sensory modalities, offering potential applications in understanding electroreception in bees, and how plants might use flowers to exploit this sensing modality. By focusing on the mathematical structure of these interactions, this work provides insight into the mechanisms underlying plant-insect communication.
For many hibernating species, the regulation of body functions is thus heavily influenced by climate cues such as temperature and snow cover. Due to the difference between adult and juvenile physiology, global warming can have qualitatively different impacts on adult and juvenile hibernator survival. To investigate the effect of climate change on hibernators' population dynamics, we consider a consumer-resource system with a continuously breeding resource and seasonally breeding hibernating consumers. In a simplistic setting, our model captures seasonality (summer and winter) with the population structured into non-reproductive juveniles and reproductive adults, while remaining analytically tractable. Our result suggests that species with faster life histories are more strongly affected by global warming. In addition, assumptions on the qualitative shape of the stage-dependent effect of global warming on hibernator winter survival can lead to qualitatively different behaviour in long-term system dynamics.
As temperatures rise, photosymbiotic marine species are presented with unique challenges. Exaptasia diaphana is a highly adaptable model organism for studying these challenges. Current methods for acquiring experimental data are expensive and require sacrificing the organism. This work focuses on the development of computational tools that will reduce cost and automate the acquisition of data. A pipeline consisting of a convolutional neural network and a transformer-based model is used to accomplish this task. Given input images of aiptasia colonies, this pipeline automatically produces accurate segmentations of aiptasia that can be used for experimental data acquisition (e.g. obtaining counts, measuring oral disk size and color information, etc.).
In an information-processing investment game, such as the growth of a population of organisms in a changing environment, Kelly betting maximizes the expected log rate of growth. In this work, we show that Kelly bets are closely related to optimal single-letter codes (i.e., they can achieve the rate-distortion bound with equality). Thus, natural information processing systems with limited computational resources can achieve information-theoretically optimal performance. We show that the rate-distortion tradeoff for an investment game has a simple linear bound, and that the bound is achievable at the point where the corresponding single-letter code is optimal. This interpretation has two interesting consequences. First, we show that increasing the organism's portfolio of potential strategies can lead to optimal performance over a continuous range of channels, even if the strategy portfolio is fixed. Second, we show that increasing an organism's number of phenotypes (i.e., its number of possible behaviours in response to the environment) can lead to higher growth rate, and we give conditions under which this occurs. Examples illustrating the results in simplified biological scenarios are presented.
Bivalve Transmissible Neoplasia (BTN) is an increasingly prevalent cancer spreading among bivalve species worldwide. BTN dynamics introduce complexities not common in many other infectious diseases due to its marine environment. Therefore, little is known about its transmission dynamics and population effects. This study develops a Bayesian framework to model BTN spread within Mya arenaria populations to address key gaps in our understanding, with a focus on statistical methodology for parameter inference and model development. We use a Bayesian compartmental modeling approach to infer and refine model parameters and leverage controlled laboratory and survey data. Laboratory data provide information on cancer cell emission, disease progression and environmental factor effects on the disease. Survey data from the field includes samples from East Coast sites that are used for fitting the model to disease progression over time in its natural environment. To capture the intricacies of BTN, our model framework builds on the traditional Susceptible-Exposed-Infectious (SEI) epidemiological model by incorporating cancer particle survivability components and the environmental effects of temperature. The Bayesian approach in our model development, implemented in STAN, provides parameter inferences and quantifies uncertainty in the results amidst limited or noisy ecological data. Future work will validate the model by comparing its predictions with 2025 survey data and conducting sensitivity analyses to identify key parameters. This statistical framework not only advances our understanding of BTN, but also demonstrates the applicability of Bayesian modeling in developing complex ecological and epidemiological compartmental models.
Wetlands are characterised by the interaction of soil with seasonal or permanent water bodies and serve several crucial ecological functions including flood prevention and water filtration. Moreover, serving as a boundary between land and aquatic environments, they provide diverse ecosystems for a variety of plants, animals and microbes. However, the ability of wetlands to sequester and store carbon also secures their place as the largest natural source of methane emissions. These emissions are strongly dependent on both the relative depth of the water table to the soil and the soil temperature, with submerged warm soil providing favourable anaerobic conditions for methanogenesis by microbes and detrimental to methane consumption through oxidation. Hence, wetland methane emissions are strongly susceptible to climate change, particularly changes in rainfall and temperature. We present a mathematical model to describe the stochastic movement of the water table coupled to a simple set of ODEs describing methane production, oxidation, and emission, parameterised by this water table depth. We employ this model to predict how the inherent variations in water table depth due to the soil profile leads to changing emission profiles across individual wetlands. Moreover, by exploring the sensitivity of these emissions to rainfall and temperature changes, we demonstrate during wetland conservation efforts the need to consider how climate change will influence emissions.
In this paper, a mathematical model is proposed to explain the interaction between Eastern Hemlock Trees and the invasive species Hemlock Woolly Adelgid. A system of reaction diffusion equations is used for this modeling exercise. There are at most three (3) steady states for the system of which the coexistent state is the only stable steady states for some parameter values. The model dynamics show that the solutions exhibit traveling wave solutions. In addition, a sensitivity analysis is conducted to determine the impact of model parameters. Sensitivity analysis suggests that the mortality rate of Eastern Hemlock Trees and the predation intensity of Hemlock Woolly Adelgid drive the dynamics of the interaction. Since Eastern Hemlock trees are foundation tress that provide shelter for several species they need to be protected for as long as possible. Based on the model dynamics and sensitivity analysis, it is postulated that selective and strategic removal of these trees will help curtail their destruction.
In the present work, using the theory of dynamical system, we discuss the dynamics of a cannibalistic predator prey model in the presence of linear harvesting schemes in prey (pest) population and provision of additional food to predators (natural enemies). A detailed mathematical analysis and numerical evaluations have been presented to discuss the pest free state, coexistence of species, stability, occurrence of different bifurcations (saddle-node, transcritical, Hopf, Bogdanov-Takens) and the impact of additional food and harvesting schemes on the dynamics of the system. Interestingly, we also observe that the pest population density decreases immediately even when small amount of harvesting is implemented. Also the eradication of pest population (stable pest free state) could be achieved via variation in the additional food and implemented harvesting schemes. The individual effects of harvesting parameters on the pest density suggest that the linear harvesting scheme is more effective to control the pest population rather than constant and nonlinear harvesting schemes.
Our work focuses on modelling how predator-prey dynamics are affected by the frequency and impact of forest disturbances such as forestry and wildfires and the subsequent regeneration of the forest. Our focus is on the Goshawk- Pine Squirrel predator-prey system. The coastal subspecies of the Goshawk is listed as endangered and their preferred habitat, mature, dense coniferous forests, is commonly disturbed. Climate change is increasing the frequency and severity of wildfires as well as periods of drought and heavy rainfall, leading to new challenges for sustaining Goshawk habitat. We present a modelling approach examining the dynamics of Goshawks and Squirrels within a disturbed landscape, and determine how different types of disturbance affect persistence of both species. An additional complicating factor that we take into account is the territoriality of Squirrels and Goshawks. Ultimately, we hope that the final model will inform conservation efforts for the Goshawk and, more generally, for predator-prey systems that rely on mature forest.
CD8+ T cells, also known as cytotoxic T cells, play a crucial role in fighting cancer by directly targeting and eliminating tumour cells. However, prolonged exposure to tumour antigens drives these cells into exhaustion, leading to the loss of their cytotoxic functions and subsequent tumour progression. The differentiation pathway undertaken by CD8+ T cells significantly influences the efficacy and persistence of the anti-tumour response. This pathway is shaped by collective inter- and intracellular decision-making processes within a complex dynamic network, involving interactions among various immune cell populations through direct cell-cell contact or signalling molecules such as cytokines. A mechanistic understanding of CD8+ T cell differentiation into specific phenotypic subsets, as well as the complex network governing this process, is essential. To address this, we develop a quantitative, data-driven mathematical model of CD8+ T cell population dynamics in response to cancer cells, capturing cell-cell interactions, cell proliferation, and T cell differentiation into effector or exhausted subsets. We analyse multiple possible network motifs governing CD8+ T cell differentiation and proliferation. In addition, we incorporate a response-time modelling approach, where the waiting-time distribution between cell states is described by a gamma rather than an exponential distribution. This approach accounts for the system’s intracellular networks in an input-to-output formulation while keeping the model’s complexity relatively manageable for analysis.
Vaccinations have historically proven to be an effective means of conferring immunity in case of various diseases by enhancing the body’s preparedness for future infection events. The success of a vaccination program depends on various factors like dose composition and time gap between vaccinations. To produce an effective response, the immune system relies heavily on B cells, among other immune cells, as these cells mature to produce antibodies. In this presentation I will present a simplified mechanistic model of B cell evolution (mutation and selection) during the immune response to vaccination, which explicitly includes the germinal centre and extrafollicular pathways. We apply our model to build an understanding of how these pathways might work together to generate a signature in the evolutionary history of B cell clonal families within a single person, considering different possible vaccination systems (homologous and heterologous). We also plan on comparing phylogenetic trees generated by our model with real trees obtained from longitudinal studies.
Synovial flares in palindromic rheumatism (PR) are aperiodic inflammatory episodes occurring in the joints, that are thought to follow a relapsing-remitting pattern. The transient and unpredictable nature of such flares is consistent with asymptomatic and non-periodic intervals. We examine the cytokine dynamics in a two-dimensional model of rheumatoid arthritis (RA) and characterise such flares as an excitable trajectory, arising out of stochastic triggers. We address questions pertaining to the frequency, decay, and persistence of synovial flares in individuals with palindromic disease. Our findings demonstrate how adaptive regulations can rescue flares that become “locked” in a 'metastable' state. However, if repetitive locking events occur over a longer timescale, they can activate a secondary adaptation toward a healthy state, which may eventually become maladaptive. Therefore, we argue that the primary mechanism underlying the progression to chronicity lies in the conflict between adaptation and maladaptation, which drives the system toward the fully developed state of rheumatoid arthritis.
BCG vaccine is the only licensed vaccine against tuberculosis (TB), a disease caused by Mycobacterium tuberculosis (Mtb). Even though billions of individuals have been vaccinated with BCG, efficacy of BCG vaccine and mechanisms by which it provides protection remain poorly understood. In a recent study, Plumlee et al. (Plos Pathogens, 19(11), e1011825, 2023) infected over a thousand mice, about half of which were vaccinated with BCG, with an ultra-low dose of Mtb (about 1 bacterium/mouse). Motivated from their study, we developed several alternative mathematical models describing Mtb dynamics in the initially infected lung (named Lung 1) and Mtb dissemination to the collateral lung (Lung 2) and fitted these models to the data from Plumlee et al. Experiments. Interestingly, proposed alternative models assuming direct or indirect Mtb dissemination describe the data well on Mtb dynamics in unvaccinated mice with similar quality. Further, we predict that Mtb replicates rapidly early during the infection, is controlled 1-2 months post-infection, and resumes replication in the chronic phase. By fitting the models to Mtb dissemination data in BCG-vaccinated mice we found that the data are best explained if BCG reduces both the rate of Mtb replication in the lungs (by 9%) and the rate of Mtb dissemination between the lungs (by 89%). Moreover, we implemented stochastic simulations of Mtb dissemination in unvaccinated and BCG-vaccinated mice, but these simulations did not fully account for the observed variability. However, stochastically simulating Mtb infection of right and left lung and dissemination between the lungs over time could successfully explain large CFU variability. Further, power analysis predicts the number of mice required in each mice group to obtain 80% power with different vaccine efficiencies. So, our mathematical modelling approach can be used to rigorously quantify efficacy of other TB vaccines in settings of ultra-low dose Mtb infection.
Correlates of protection against infection with Mycobacterium tuberculosis (Mtb) or against tuberculosis (TB) remain poorly defined. The ratio of colony forming units (CFUs) to chromosomal equivalents (CEQs), Z = CFU/CEQ, has been used as a metric for how effectively Mtb is killed in vivo. However, the contribution of bacterial killing to changes in CFU/CEQ ratio during an infection has not been rigorously investigated. We developed alternative mathe- matical models to study the dynamics of CFUs, CEQs, and Z during an Mtb infection. We find that the ratio Z alone cannot be used to infer the death rate of bacteria, unless the dynamics of CEQs and CFUs are entirely uncoupled, which is biologically unreasonable and inconsistent with the view that CEQs reflect an accumulation of both viable and non-viable bacteria. We estimate a half life of about 20 days of Mtb H37Rv CEQs in mice, similar to that found for Mtb Erdman in cynomolgus macaques. Although this seems slow, we found that estimated rates of Mtb replication and death are extremely sensitive even to slow decay of detectable Mtb genomes. We provide evidence of substantial killing of Mtb bacteria prior to arrival of adaptive immunity to the site of infection. We also propose experiments that will allow to more accurately measure the rate of Mtb DNA loss helping more rigorously to quantify impact of immunity on within-host Mtb dynamics.
Chronic inflammation disrupts hematopoietic homeostasis, causing pathological myelopoiesis and malignant clones that grow. The study uses a mathematical model with local (bone marrow) and global (peripheral inflammation) negative feedback mechanisms to examine how inflammation-driven regulations affect HSC self-renewal, progenitor dynamics, and differentiation. Healthy and malignant populations compete in the model, which examines system stability through feedback mechanisms. The results show that chronic inflammation can cause myelopoietic disorders by overproducing progenitor cells and disrupting lineage balance without global feedback regulation. Self-renewal feedback regulates stem cell proliferation to strengthen hematopoietic cells and mitigate chronic inflammation damages. Because excessive suppression can destabilize hematopoiesis, the model suggests tightly controlling negative feedback on progenitor cells. Mutations affecting global feedback can cause malignant clones, revealing how inflammation causes hematological malignancies like MDS and AML.
Part of the immune response upon infection involves B cells and a process known as affinity maturation. During affinity maturation, produced antibodies increase in affinity to presented antigen. Additionally, plasma B cells and memory B cells are created. This is to allow the system to remember and quickly mount a response to the presented antigen in the case of a repeat infection. Repeated exposures to the same antigen will produce antibodies of successively greater affinities. However, as antigen move away in antigenic distance from the initial strain (antigenic drift), the ability of the body to cross-reactively neutralize the antigen decreases. This issue has been well documented in cases of influenza and there is a concern it is occurring in SARS-CoV-2 given successive variants of concern (VOC). Such VOCs would be less susceptible to any immune protection gained from vaccination and prior infection. We modeled these processes using an agent-based model (ABM) that considers B cells (naïve, plasma, memory), antibodies, and antigens. We represent receptor (B cells, antibodies) and epitope (antigens) proteins in Euclidean shape space, simulating binding between these agents based on Hamming distance. We also consider the formation of immune complexes—free antibodies bound to antigen which limits the antigen’s ability to infect more cells. We simulated SARS-CoV-2 infections using our ABM. We present results that examine immune responses when presented with various VOCs and differing immune imprinting.
Some viruses exhibit “rebound” when the administration of antiviral drugs is discontinued. Viral rebound caused by resistance mutations or latent reservoirs has been studied mathematically. In this study, we investigated the viral rebound due to other causes. Since immunity is weaker during antiviral treatment than without the treatment, drug discontinuation may lead to an increase in the viral load. We analyzed the dynamics of the number of virus-infected cells, cytotoxic T lymphocytes, and memory cells and identified the conditions under which the viral load increased upon drug discontinuation. If drug is administered for an extended period, a viral rebound occurs when the ratio of viral growth rate in the absence to that in the presence of the antiviral drug exceeds the “rebound threshold.” We analyzed how the rebound threshold depended on the patient’s conditions and the type of treatment. Mathematical and numerical analyses revealed that rebound after discontinuation was more likely to occur when the drug effectively reduced viral proliferation, drug discontinuation was delayed, and the processes activating immune responses directly were stronger than those occurring indirectly through immune memory formation. We discussed additional reasons for drugs to cause viral rebound more likely.
COVID-19 outcomes vary widely among individuals, with most having mild illness, while a small percentage experience severe symptoms and a minor fraction death. Several treatments for COVID-19 have been proposed. One of the most promising is the inhibition of mTORC1 by Sirolimus. However, not all patients are sensitive to this treatment. To uncover the complex relations behind the heterogeneity and sensitivity of some individuals to treatments, we developed a hybrid agent-based model of the innate immune response to study the infection in the whole lung. The model includes key cells involved in the disease and critical intracellular factors such as NF-kB, IRF3, STAT1, and mTORC1. We calibrated and validated our model using literature and our own experimental data. We used it to explore different scenarios and explain our experimental results showing a positive correlation between mTORC1 activity and viral replication but a negative correlation between mTORC1 and IFN-a expression. Our initial simulations showed that mTORC1 is a master regulator of intracellular viral response and suggested novel intervention targets upstream of mTORC1. Our aim is to personalize the model and quantify the role of mTORC1 in the COVID-19 heterogeneity.
CD8+ T cells, also known as cytotoxic T cells, play a crucial role in fighting cancer by directly targeting and eliminating tumour cells. However, prolonged exposure to tumour antigens drives these cells into exhaustion, leading to the loss of their cytotoxic functions and subsequent tumour progression. The differentiation pathway undertaken by CD8+ T cells significantly influences the efficacy and persistence of the anti-tumour response. This pathway is shaped by collective inter- and intracellular decision-making processes within a complex dynamic network, involving interactions among various immune cell populations through direct cell-cell contact or signalling molecules such as cytokines. A mechanistic understanding of CD8+ T cell differentiation into specific phenotypic subsets, as well as the complex network governing this process, is essential. To address this, we develop a quantitative, data-driven mathematical model of CD8+ T cell population dynamics in response to cancer cells, capturing cell-cell interactions, cell proliferation, and T cell differentiation into effector or exhausted subsets. We analyse multiple possible network motifs governing CD8+ T cell differentiation and proliferation. In addition, we incorporate a response-time modelling approach, where the waiting-time distribution between cell states is described by a gamma rather than an exponential distribution. This approach accounts for the system’s intracellular networks in an input-to-output formulation while keeping the model’s complexity relatively manageable for analysis.
In this paper, we study the immune system’s response to infection with the bacteria Mycobacterium tuberculosis (the causative agent of tuberculosis). The response by the immune system is either global (lymph node, thymus, and blood) or local (at the site of infection). The response by the immune system against tuberculosis (TB) at the site of infection leads to the formation of spherical structures which comprised of cells, bacteria, and effector molecules known as granuloma. We developed a deterministic model capturing the dynamics of the immune system, macrophages, cytokines and bacteria. The hallmark of Mycobacterium tuberculosis (MTB) infection in the early stages requires a strong protective cell-mediated naive T cells differentiation which is characterised by antigen-specific interferon gamma (IFN-γ). The host immune response is believed to be regulated by the interleukin-10 cytokine by playing the critical role of orchestrating the T helper 1 and T helper 2 dominance during disease progression. The basic reproduction number is computed and a stability analysis of the equilibrium points is also performed. Through the computation of the reproduction number, we predict disease progression scenario including the latency state. The occurrence of latent infection is shown to depend on a number of effector function and the bacterial load for R0 < 1. The model predicts that endemically there is no steady state behaviour; rather it depicts the existence of the MTB to be a continuous process progressing over a differing time period. Simulations of the model predict the time at which the activated macrophages overcome the infected macrophages (switching time) and observed that the activation rate (ω) correlates negatively with it. The efficacy of potential host-directed therapies was determined by the use of the model.
Capsular contracture is a pathological response to implant-based reconstructive breast surgery, where the ``capsule'' (tissue surrounding an implant) painfully thickens, contracts and deforms. It is known to affect breast-cancer survivors at higher rates than healthy women opting for cosmetic breast augmentation with implants. We model the early stages of capsular contracture based on stress-dependent recruitment of contractile and mechanosensitive cells to the implant site. We derive a one-dimensional continuum spatial model for the spatio-temporal evolution of cells and collagen densities away from the implant surface. Various mechanistic assumptions are investigated for linear versus saturating mechanical cell responses and cell traction forces. Our results point to specific risk factors for capsular contracture, and indicate how physiological parameters, as well as initial states (such as inflammation after surgery) contribute to patient susceptibility.
Capsular contracture is a pathological response to implant-based reconstructive breast surgery, where the ``capsule'' (tissue surrounding an implant) painfully thickens, contracts and deforms. It is known to affect breast-cancer survivors at higher rates than healthy women opting for cosmetic breast augmentation with implants. We model the early stages of capsular contracture based on stress-dependent recruitment of contractile and mechanosensitive cells to the implant site. We derive a one-dimensional continuum spatial model for the spatio-temporal evolution of cells and collagen densities away from the implant surface. Various mechanistic assumptions are investigated for linear versus saturating mechanical cell responses and cell traction forces. Our results point to specific risk factors for capsular contracture, and indicate how physiological parameters, as well as initial states (such as inflammation after surgery) contribute to patient susceptibility.
Neuroinflammation is driven by cellular interactions that are difficult to capture with experimental or mathematical approaches alone. Microglia, the resident immune cells of the central nervous system, dynamically shift between functional states in response to different stimuli. Here we present an integrated framework that combines microscopy, cell quantification, and RNA sequencing with mathematical modeling to describe the cellular interactions underlying neuroinflammation. Our data reveal four distinct microglia subtypes whose behaviors and interactions with other neural cell types are incorporated into a system of ODEs. This approach allows us to move from population-level description of microglia to a subtype-specific description of their roles in inflammation. Using our mathematical model, we explore how microglia subpopulations differentially contribute to the propagation of the neuroimmune response, specifically through cytokine secretion, blood-brain barrier weakening, and T cell activation. We will also outline how our suite of primary experimental data can be used by others to inform future mathematical models.
Here I report on work with Yogesh Bali on response functions in models for T cell activation. How does the activity of T cells depend on the amount of antigen they are exposed to and the dissociation time of the binding of the antigen to the T cell receptor? We have studied a number of ODE models addressing this question using analytical and numerical techniques. Here I highlight two key results of this work. First, it can happen that for biologically reasonable parameter values an increase in the dissociation time can lead to a decrease in the response. Second, it can happen that the dependence of response on the control parameters exhibits more than one maximum.
Invasive pulmonary aspergillosis is a deadly disease caused by the mold Aspergillus. As the mold grows in the lungs, fungal hyphae penetrate the epithelium, resulting in lung hemorrhage. We previously reported that extracellular heme worsens the outcome of the infection. We hypothesize that hemostasis is protective in invasive aspergillosis. Methods: C57Bl/6J mice were partially neutrophil-depleted and challenged with Aspergillus conidia. We performed serial thromboelastography on the blood of infected mice and control mice, sampled the alveolar lumen by bronchoalveolar lavage (BAL. We also performed ELISAs for coagulation factor Xa and Thrombin-antithrombin complex on BAL. We used mathematical modeling to map coagulation factors to thromboelastography curves and predict coagulation factors that explain observed TEG patterns. In some experiments, infected mice were treated with clinical drugs that inhibit factor Xa (apixaban) and prevent fibrinolysis (tranexamic acid). Results: Infected mice had higher levels of factor Xa and thrombin-antithrombin complex in BAL, and higher maximum amplitudes in thromboelastography compared to uninfected mice, indicating appropriate activation of the coagulation. Unexpectedly, infected animals had an elongated time to clot on thromboelastography. Our model predicted a potential depletion of factors X or VII. We found a partial depletion of factor VII but not factor X in the blood. Treatment with apixaban increased fungal burden. Treatment with tranexamic acid resulted in a pronounced reduction in fungal burden in female mice but had no effect on male mice. Conclusions: Our preliminary studies suggest that coagulopathy is an important component of invasive aspergillosis and that treatment with anticoagulants during infection might lead to worse outcomes in mice. Treatment with an anti-fibrinolytic agent led to a lowered fungal burden in female mice, and agents that aid in clot formation might improve outcomes in mice.
Understanding how people living with HIV (PLWH) respond to repeated COVID-19 vaccinations is critical for advancing precision medicine in immunocompromised populations. In this study, we use random forest models to identify which immune responses most effectively differentiate vaccine outcomes between PLWH on antiretroviral therapy and an HIV-negative control group. Our data set contains an extensive range of immune features, including serum and saliva IgG and IgA responses, ELISpot IFNg and IL2 responses to SARS-CoV-2 spike peptides, ACE2 receptor displacement, and SARS-CoV-2 neutralization capacity; all tracked longitudinally up to 104 weeks in each individual following SARS-CoV-2 vaccine dose 1, up to dose 5. Our model achieves near-perfect accuracy and reveals that cytokine-producing T cells and saliva-based IgA responses are key features for classification, whereas serum IgG markers offer limited classification value. Through ablation sensitivity analysis, we are able to identify outlier HIV- and HIV+ individuals whose immunological profiles do not fit the learned ‘pattern’ identified by the RF algorithm; some HIV+ individuals on ART appear to have nearly complete immune recovery while some HIV- individuals have vaccine-elicited immune signatures that appear like that of a typical HIV+ individual, suggesting immune compromisation.
NF-κB signaling shapes the inflammatory response, and its dysregulation is linked with autoimmune, neurodegenerative, and cardiovascular diseases. After pathway activation, dimers of the NF-κB family act as transcription factors for a large set of target genes including positive and negative regulators of the upstream pathway. Important examples are cytokines, such as TNFα, that stimulates the pathway and therefore contributes to the intercellular communication of cells. While the cellular NF-κB pathway has been intensively studied via computational modeling the effect of intercellular coupling is less explored. To investigate this in spatial contexts, we build a multi-scale, multi-cell ordinary differential equation model where physiological cell properties are computed with a Cellular Potts Model using the Morpheus software. To elucidate the interaction between different cell types in the intestinal crypt we extend the model to two NF-κB expressing cell types: sentinel macrophages and epithelial cells. There, the LPS-activated macrophages secrete TNFα to elicit an immune response in the epithel. We evaluated the models in terms of sensitivity regarding signal transmission strength and speed as well as common single-cell measures, like maximum NF-κB activation and time thereoff.
Vaccinations have historically proven to be an effective means of conferring immunity in case of various diseases by enhancing the body’s preparedness for future infection events. The success of a vaccination program depends on various factors like dose composition and time gap between vaccinations. To produce an effective response, the immune system relies heavily on B cells, among other immune cells, as these cells mature to produce antibodies. In this presentation I will present a simplified mechanistic model of B cell evolution (mutation and selection) during the immune response to vaccination, which explicitly includes the germinal centre and extrafollicular pathways. We apply our model to build an understanding of how these pathways might work together to generate a signature in the evolutionary history of B cell clonal families within a single person, considering different possible vaccination systems (homologous and heterologous). We also plan on comparing phylogenetic trees generated by our model with real trees obtained from longitudinal studies.
Synovial flares in palindromic rheumatism (PR) are aperiodic inflammatory episodes occurring in the joints, that are thought to follow a relapsing-remitting pattern. The transient and unpredictable nature of such flares is consistent with asymptomatic and non-periodic intervals. We examine the cytokine dynamics in a two-dimensional model of rheumatoid arthritis (RA) and characterise such flares as an excitable trajectory, arising out of stochastic triggers. We address questions pertaining to the frequency, decay, and persistence of synovial flares in individuals with palindromic disease. Our findings demonstrate how adaptive regulations can rescue flares that become “locked” in a 'metastable' state. However, if repetitive locking events occur over a longer timescale, they can activate a secondary adaptation toward a healthy state, which may eventually become maladaptive. Therefore, we argue that the primary mechanism underlying the progression to chronicity lies in the conflict between adaptation and maladaptation, which drives the system toward the fully developed state of rheumatoid arthritis.
BCG vaccine is the only licensed vaccine against tuberculosis (TB), a disease caused by Mycobacterium tuberculosis (Mtb). Even though billions of individuals have been vaccinated with BCG, efficacy of BCG vaccine and mechanisms by which it provides protection remain poorly understood. In a recent study, Plumlee et al. (Plos Pathogens, 19(11), e1011825, 2023) infected over a thousand mice, about half of which were vaccinated with BCG, with an ultra-low dose of Mtb (about 1 bacterium/mouse). Motivated from their study, we developed several alternative mathematical models describing Mtb dynamics in the initially infected lung (named Lung 1) and Mtb dissemination to the collateral lung (Lung 2) and fitted these models to the data from Plumlee et al. Experiments. Interestingly, proposed alternative models assuming direct or indirect Mtb dissemination describe the data well on Mtb dynamics in unvaccinated mice with similar quality. Further, we predict that Mtb replicates rapidly early during the infection, is controlled 1-2 months post-infection, and resumes replication in the chronic phase. By fitting the models to Mtb dissemination data in BCG-vaccinated mice we found that the data are best explained if BCG reduces both the rate of Mtb replication in the lungs (by 9%) and the rate of Mtb dissemination between the lungs (by 89%). Moreover, we implemented stochastic simulations of Mtb dissemination in unvaccinated and BCG-vaccinated mice, but these simulations did not fully account for the observed variability. However, stochastically simulating Mtb infection of right and left lung and dissemination between the lungs over time could successfully explain large CFU variability. Further, power analysis predicts the number of mice required in each mice group to obtain 80% power with different vaccine efficiencies. So, our mathematical modelling approach can be used to rigorously quantify efficacy of other TB vaccines in settings of ultra-low dose Mtb infection.
Correlates of protection against infection with Mycobacterium tuberculosis (Mtb) or against tuberculosis (TB) remain poorly defined. The ratio of colony forming units (CFUs) to chromosomal equivalents (CEQs), Z = CFU/CEQ, has been used as a metric for how effectively Mtb is killed in vivo. However, the contribution of bacterial killing to changes in CFU/CEQ ratio during an infection has not been rigorously investigated. We developed alternative mathe- matical models to study the dynamics of CFUs, CEQs, and Z during an Mtb infection. We find that the ratio Z alone cannot be used to infer the death rate of bacteria, unless the dynamics of CEQs and CFUs are entirely uncoupled, which is biologically unreasonable and inconsistent with the view that CEQs reflect an accumulation of both viable and non-viable bacteria. We estimate a half life of about 20 days of Mtb H37Rv CEQs in mice, similar to that found for Mtb Erdman in cynomolgus macaques. Although this seems slow, we found that estimated rates of Mtb replication and death are extremely sensitive even to slow decay of detectable Mtb genomes. We provide evidence of substantial killing of Mtb bacteria prior to arrival of adaptive immunity to the site of infection. We also propose experiments that will allow to more accurately measure the rate of Mtb DNA loss helping more rigorously to quantify impact of immunity on within-host Mtb dynamics.
Chronic inflammation disrupts hematopoietic homeostasis, causing pathological myelopoiesis and malignant clones that grow. The study uses a mathematical model with local (bone marrow) and global (peripheral inflammation) negative feedback mechanisms to examine how inflammation-driven regulations affect HSC self-renewal, progenitor dynamics, and differentiation. Healthy and malignant populations compete in the model, which examines system stability through feedback mechanisms. The results show that chronic inflammation can cause myelopoietic disorders by overproducing progenitor cells and disrupting lineage balance without global feedback regulation. Self-renewal feedback regulates stem cell proliferation to strengthen hematopoietic cells and mitigate chronic inflammation damages. Because excessive suppression can destabilize hematopoiesis, the model suggests tightly controlling negative feedback on progenitor cells. Mutations affecting global feedback can cause malignant clones, revealing how inflammation causes hematological malignancies like MDS and AML.
Part of the immune response upon infection involves B cells and a process known as affinity maturation. During affinity maturation, produced antibodies increase in affinity to presented antigen. Additionally, plasma B cells and memory B cells are created. This is to allow the system to remember and quickly mount a response to the presented antigen in the case of a repeat infection. Repeated exposures to the same antigen will produce antibodies of successively greater affinities. However, as antigen move away in antigenic distance from the initial strain (antigenic drift), the ability of the body to cross-reactively neutralize the antigen decreases. This issue has been well documented in cases of influenza and there is a concern it is occurring in SARS-CoV-2 given successive variants of concern (VOC). Such VOCs would be less susceptible to any immune protection gained from vaccination and prior infection. We modeled these processes using an agent-based model (ABM) that considers B cells (naïve, plasma, memory), antibodies, and antigens. We represent receptor (B cells, antibodies) and epitope (antigens) proteins in Euclidean shape space, simulating binding between these agents based on Hamming distance. We also consider the formation of immune complexes—free antibodies bound to antigen which limits the antigen’s ability to infect more cells. We simulated SARS-CoV-2 infections using our ABM. We present results that examine immune responses when presented with various VOCs and differing immune imprinting.
Some viruses exhibit “rebound” when the administration of antiviral drugs is discontinued. Viral rebound caused by resistance mutations or latent reservoirs has been studied mathematically. In this study, we investigated the viral rebound due to other causes. Since immunity is weaker during antiviral treatment than without the treatment, drug discontinuation may lead to an increase in the viral load. We analyzed the dynamics of the number of virus-infected cells, cytotoxic T lymphocytes, and memory cells and identified the conditions under which the viral load increased upon drug discontinuation. If drug is administered for an extended period, a viral rebound occurs when the ratio of viral growth rate in the absence to that in the presence of the antiviral drug exceeds the “rebound threshold.” We analyzed how the rebound threshold depended on the patient’s conditions and the type of treatment. Mathematical and numerical analyses revealed that rebound after discontinuation was more likely to occur when the drug effectively reduced viral proliferation, drug discontinuation was delayed, and the processes activating immune responses directly were stronger than those occurring indirectly through immune memory formation. We discussed additional reasons for drugs to cause viral rebound more likely.
COVID-19 outcomes vary widely among individuals, with most having mild illness, while a small percentage experience severe symptoms and a minor fraction death. Several treatments for COVID-19 have been proposed. One of the most promising is the inhibition of mTORC1 by Sirolimus. However, not all patients are sensitive to this treatment. To uncover the complex relations behind the heterogeneity and sensitivity of some individuals to treatments, we developed a hybrid agent-based model of the innate immune response to study the infection in the whole lung. The model includes key cells involved in the disease and critical intracellular factors such as NF-kB, IRF3, STAT1, and mTORC1. We calibrated and validated our model using literature and our own experimental data. We used it to explore different scenarios and explain our experimental results showing a positive correlation between mTORC1 activity and viral replication but a negative correlation between mTORC1 and IFN-a expression. Our initial simulations showed that mTORC1 is a master regulator of intracellular viral response and suggested novel intervention targets upstream of mTORC1. Our aim is to personalize the model and quantify the role of mTORC1 in the COVID-19 heterogeneity.
A metapopulation consists of a group of spatially distanced subpopulations, each occupying a separate patch. It is usually assumed that each localized patch is well-mixed. In this talk, we will discuss the spread of an epidemic in a system of weakly connected patches, where the disease dynamics of each patch occurs on a network. The SIR dynamics in a single patch is governed by the rate of disease transmission, the disease duration, and the node degree distribution of a network. Monte-Carlo simulations of the model reveal the phenomenon of spatial disease propagation. The speed of front propagation and its dependence on the single patch parameters and on the strength of interaction between the patches was determined analytically, and a good agreement with simulation results was observed [1]. Next, we will discuss front propagation in case of an Allee effect, where the effective transmission rate depends on the fraction of infected, and the state of no epidemic is linearly stable. We discovered [2] a novel phenomenon of front stoppage: in some regime of parameters, the front solution ceases to exist, and the propagating pulse of infection decays despite the initial outbreak. [1]. E. Khain and M. Iyengar, Phys. Rev. E 107, 034309 (2023). [2]. E. Khain, Phys. Rev. E 107, 064303 (2023).
Stochasticity is introduced to bring new insight into the modelling of population dynamics of diseases. Many systems, in nature, are subject to stochastic perturbations. In this talk, we present differential equations with stochastic perturbations and the updated data estimation method for estimating the transmission rate changes over time. The models for the population dynamics under SEIR epidemic models with stochastic perturbations are presented for the dynamics of the COVID-19 pandemic in Bogota, Colombia. We performed computational experiments to interpret COVID-19 dynamics using actual data for the proposed models. We estimated the model parameters and updated their reported infected and recovered data estimates. (joint work with Andres Rios-Gutierrez )
Self-isolation and stay-at-home measures are crucial for curbing the spread of contagious pathogens while vaccines are being developed. Furthermore, research during the 2019-22 coronavirus pandemic (COVID-19) emphasizes that proper enforcement and timely lifting of these measures are vital for effective disease management. In this context, we analyzed a simple dynamical system to understand how an epidemic progresses by isolating susceptible individuals (confinement) and reintroducing them to infection (deconfinement). This model captures the overall magnitude and direction of flows between confined and deconfined groups—akin to osmosis—leading to a dimensionless quantity defined as confinement tonicity. Our mathematical analysis suggests that confinement tonicity influences the final epidemic size, providing insights into careful quarantine management for effective disease control.
The accrual of nucleotide substitutions in pathogen genomes accompanies their transmission through host populations. Because lineages with higher fitness tend to transmit rapidly to new hosts before incurring very many substitutions, large numbers of related sequences are usually interpreted as evidence of transmission success. Quantities like the local branching index (LBI) aim to identify successful lineages in this way by scoring sequences according to the number of close relatives captured in the dataset. While statistics like LBI are easily calculated from a given phylogenetic tree (or a distribution of trees), observation errors related to sampling bias and censoring may introduce spurious signals of transmission success. To disentangle these effects, we use stochastic compartmental models to simulate outbreaks and generate distributions of phylogenies under a variety of testing programs (such as surveillance of symptomatic cases, or cross-sectional prevalence studies). By characterizing the types of phylogenies expected under these situations, we can work towards a clearer understanding of the types of signals that are likely to be detected with sequence data.
Dengue fever has emerged as an increasingly alarming public health challenge, further complicated by the impacts of climate change on control efforts. Yet, the full extent of climate's impact on dengue incidence remains poorly understood. To investigate this, we employed an advanced causal inference method to 16 regions in the Philippines, selected for their diverse climatic conditions. Unlike previous methods for detecting regulatory relationships, this method is capable of detecting nonlinear and joint effects of temperature and rainfall to dengue incidence. We found that temperature consistently increased dengue incidence throughout all the regions, while rainfall effects differed depending on location. Further analysis showed that this pattern is due to the variation in dry season length, a factor previously overlooked. Specifically, our results showed that regions with low variation in dry season length experience a negative impact of rainfall on dengue incidence likely due to strong flushing effect on mosquito habitats, while regions with high variation in dry season length experience a positive impact, likely due to increased mosquito breeding sites. This study offers a fresh perspective on the relationship between climate and dengue incidence, emphasizing the need for tailored prevention strategies based on local climate conditions.
This study investigates the spatial-temporal heterogeneity in the relationship between human mobility and COVID-19 transmission across 229 regions in South Korea during six epidemic waves from January 2020 to September 2022. While previous research primarily focused on the early stages of the pandemic and the impacts of mobility restrictions, our study utilizes mobility data from SK Telecom and COVID-19 case data from the Korea Disease Control and Prevention Agency to provide a more comprehensive analysis. We applied empirical mode decomposition (EMD) and clustering analysis to classify regional mobility patterns and conducted cross-correlation analysis to assess the relationship between mobility and confirmed cases. The findings indicate that incoming mobility significantly influenced the number of confirmed cases in urban and densely populated areas, whereas rural regions exhibited contrasting patterns. Moreover, these relationships evolved across different epidemic waves, highlighting the influence of regional characteristics and public health interventions. This study underscores the need to consider spatial-temporal heterogeneity in mobility-transmission dynamics to develop tailored public health strategies and enhance preparedness for future pandemics.
Forecasting the occurrence and absence of novel disease outbreaks is essential for disease management, yet existing methods are often context-specific, require a long preparation time, and non-outbreak prediction remains understudied. To address this gap, we propose a novel framework using a feature-based time series classification (TSC) method to forecast outbreaks and non-outbreaks. We tested our methods on synthetic data from a Susceptible–Infected–Recovered (SIR) model for slowly changing, noisy disease dynamics. Outbreak sequences give a transcritical bifurcation within a specified future time window, whereas non-outbreak (null bifurcation) sequences do not. We identified incipient differences, reflected in 22 statistical features and 5 early warning signal indicators, in time series of infectives leading to future outbreaks and non-outbreaks. Classifier performance, given by the area under the receiver-operating curve (AUC), ranged from 0.99 for large expanding windows of training data to 0.7 for small rolling windows. The framework is further evaluated on four empirical datasets: COVID-19 incidence data from Singapore, 18 other countries, and Edmonton, Canada, as well as SARS data from Hong Kong, with two classifiers exhibiting consistently high accuracy. Our results highlight detectable statistical features distinguishing outbreak and non-outbreak sequences well before potential occurrence, in both synthetic and real-world datasets presented in this study.
The increase in Plasmodium vivax malaria cases in Korea highlights the necessity to reevaluate intervention strategies as climate patterns change. In 2024, confirmed cases rose by 37% compared to the previous three years' average, along with an increase in vector mosquito populations. In response, the Korea Disease Control and Prevention Agency (KCDA) expanded designated malaria risk areas in Seoul. Effective control depends on optimizing non-pharmaceutical interventions with primaquine-based treatment. As tafenoquine emerges as a potential alternative treatment, evaluating its impact on malaria transmission, relapse rate and cost-effectiveness within public health systems is essential. To address these issues, we developed a mathematical model incorporating climate variability to assess the effectiveness of non-pharmaceutical interventions under different climate scenarios. Using the Improved Multi-Objective Differential Evolution (IMODE) algorithm, we analyzed the optimal interventions based on observed malaria control measures. Our results suggest that optimal intervention strategies can significantly reduce malaria transmission and relapse rate, highlighting the cost-effectiveness of tafenoquine and optimal intervention approaches in Korea’s malaria control measures.
Abstract Background: Healthcare accessibility is a key determinant of health outcomes during pandemics. Disparities in access may contribute to indirect excess mortality beyond reported COVID-19 deaths. This study quantitatively examines the impact of regional healthcare accessibility on non-COVID excess mortality in South Korea using a mathematical modeling approach. Methods: We first performed hierarchical clustering based on the average travel time to various healthcare facilities, classifying regions into two groups: Cluster 0 (high accessibility) and Cluster 1 (low accessibility). A CatBoost model trained on 2014–2019 data predicted expected deaths for 2020–2022, and excess mortality was calculated as the difference between observed and predicted deaths. Finally, multiple linear regression was then used to evaluate the association between accessibility time and non-COVID excess mortality. Results: Our analysis revealed that regions with high healthcare accessibility (Cluster 0) exhibited excess mortality patterns that closely aligned with reported COVID-19 deaths. In contrast, regions with lower accessibility (Cluster 1) experienced a significant increase in non-COVID excess mortality, particularly during the Omicron surge (fifth and sixth pandemic waves). The regression analysis demonstrated that longer healthcare accessibility times were significantly associated with higher non-COVID excess mortality in later pandemic stages. Conclusion: This study demonstrates that regional disparities in healthcare accessibility contribute to indirect excess mortality during pandemics. The findings highlight the importance of targeted policy interventions, such as strengthening healthcare infrastructure and expanding telemedicine, to reduce health inequalities and enhance public health resilience in future crises.
A metapopulation consists of a group of spatially distanced subpopulations, each occupying a separate patch. It is usually assumed that each localized patch is well-mixed. In this talk, we will discuss the spread of an epidemic in a system of weakly connected patches, where the disease dynamics of each patch occurs on a network. The SIR dynamics in a single patch is governed by the rate of disease transmission, the disease duration, and the node degree distribution of a network. Monte-Carlo simulations of the model reveal the phenomenon of spatial disease propagation. The speed of front propagation and its dependence on the single patch parameters and on the strength of interaction between the patches was determined analytically, and a good agreement with simulation results was observed [1]. Next, we will discuss front propagation in case of an Allee effect, where the effective transmission rate depends on the fraction of infected, and the state of no epidemic is linearly stable. We discovered [2] a novel phenomenon of front stoppage: in some regime of parameters, the front solution ceases to exist, and the propagating pulse of infection decays despite the initial outbreak. [1]. E. Khain and M. Iyengar, Phys. Rev. E 107, 034309 (2023). [2]. E. Khain, Phys. Rev. E 107, 064303 (2023).
Epidemic modeling is essential for understanding and managing the spread of infectious diseases. However, it often faces challenges related to unidentifiability due to high-dimensional parameters. Therefore, integrating various data sources to infer epidemic parameters is crucial for reliable modeling. We propose MPUGAT, a hybrid framework that combines a multi-patch compartmental model with a spatiotemporal deep learning approach. By leveraging a Graph Attention Network (GAT), MPUGAT effectively captures spatiotemporal infection patterns from diverse time series data to infer a dynamic transmission matrix. Applied to COVID-19 data from South Korea, MPUGAT demonstrates superior performance in estimating the time-varying transmission matrix, aligning well with real-world dynamics. This framework offers a novel approach to integrating easily accessible, low-dimensional, non-epidemic-related data into epidemic modeling, enhancing both inference and interpretability.
Dengue virus (DENV) has been causing outbreaks and epidemics over the course of the whole XX century. Recently, the size of the seasonal epidemics has been sequentially reaching record-breaking numbers, in 2023 WHO reported globally a record of over 6.8 million infections, and last year the WHO reported again a record-breaking number of confirmed cases, with over 10.5 million confirmed cases. Mainly those confirmed cases have happened in the Americas which has reportedly been affected by different serotypes of the virus. The region of Central America and the Caribbean is mainly affected by DENV3, while the South America region is affected in major number by the DENV1 and DENV2 serotypes. This situation raises the concern of how immunity and the ecological niche of the serotypes works and how this can be better understood to help design plans of contingency. To do so, we have adapted a well-known model to many-strain pathogen dynamics from the point of view of the strains. The model keeps the dynamics simple while being robust in incorporate as many as needed different strains. We modified the model to first reproduce the dengue dynamics in human and mosquitoes population, from that we change the demographics of each population to study how different relative time of infection to life span can give rises to different effects in the strain space and infection niches. All four serotypes have almost total homotypic immunity and, in the short term, heterotypic immunity. The goal here is to have a simple formulation for all the four serotypes and understand how this different dynamical regimes on different hosts affect i) the emergence of niche to the strains, and ii) how it determines endemicity of the disease in humans.
Secondary vaccine failure(SVF) following the second dose of measles-mumps-rubella(MMR) vaccine has resulted in low seroprevalence among healthcare workers(HCWs) in their 20s in the Republic of Korea. During the 2019 measles outbreak, 73% of confirmed cases in a hospital were seropositive yet still infected, highlighting that the presence of antibodies does not guarantee full protection. This study evaluates the impact of waning immunity on future measles outbreaks and develops vaccination strategies to control the nosocomial transmission. We developed a SEIR model incorporating three immunity states; Protected, Partially protected, and Seronegative, and integrated hospital seroprevalence and age structure. Using the stochastic Gillespie algorithm, we simulated the 2019 outbreak and predicted future scenarios. Our analysis revealed that the transmission rate among seronegative individuals was approximately 2.66 times higher than that of partially protected individuals. In long-term projections, vaccination only for new HCWs reduced confirmed cases by 41–51% compared to no vaccination. In contrast, vaccination for all HCWs suppressed outbreaks for approximately one year by reducing the effective reproduction number below 1. However, infections among partially protected individuals caused the overall outbreak size to increase over time. While current guidelines for third dose of MMR focus on seronegative individuals, our study provides mathematical evidence that booster shots for all HCWs are a more effective strategy than targeting only seronegative individuals in nosocomial environments.
Seasonal influenza remains a significant public health concern in South Korea, with seasonal outbreaks contributing to a high burden of disease and healthcare costs. To inform effective vaccination strategies, mathematical models can provide insights into the dynamics of influenza transmission and the impact of various vaccination strategies. In this talk, I will discuss a compartmental model to simulate the spread of seasonal influenza in South Korea, incorporating factors such as vaccination coverage, demographic structure, and virus transmission characteristics. The model is calibrated using historical data on influenza incidence, vaccination rates, and population demographics. I will discuss different vaccination strategies, including age-targeted vaccination, the impact of varying vaccine efficacy levels. The findings of considered work suggest that a targeted vaccination strategy, focusing on high-risk groups, combined with public health measures to increase vaccine uptake, could significantly reduce the incidence of influenza and mitigate its economic impact. This mathematical framework provides valuable guidance for optimizing influenza vaccination policies in South Korea and can be adapted to other countries facing similar challenges with seasonal influenza control.
Ticks are known as the important vectors that can carry and spread various tick-borne diseases by biting humans. Representative tick-borne diseases include Lyme disease, which is prevalent worldwide, tick-borne encephalitis, which is mainly prevalent in Europe and Africa, and Severe Fever Thrombocytopenia Syndrome (SFTS), one of the emerging infectious diseases spreading in East Asia. Among these tick-borne diseases, severe fever thrombocytopenia syndrome (SFTS) is an infectious disease that was first recognized in the 2010s and has recently threatened human health. Since the disease was discovered less than 20 years ago, research on the ecology of vectors and transmission dynamics of SFTS is still insufficient. Therefore, in this study, we developed the stage-structured mathematical model for the population dynamics of Haemaphysalis longicornis (Asian Longhorned Ticks), the main vector transmitting SFTS, and extended the model to include the transmission dynamics of SFTS in Korea. Based on the model, we analyzed the impact of climate change on tick population with climate-dependent parameters in the model. Additionally, the effects of control measures have been investigated following the changes in the tick population, SFTS patients, and the costs associated with SFTS.
Mpox (formerly known as monkeypox) is a neglected tropical disease that became notorious during its 2022-2023 worldwide outbreak. The vaccination was available, but there were inequities in vaccine access. In this paper, we extend existing game-theoretic models to study a population that is heterogeneous in the relative vaccination costs. We consider a population with two groups. We determine the Nash equilibria (NE), i.e., optimal vaccination rates, for each of the groups. We show that the NE always exists and that, for a narrow range of parameter values, there can be multiple NEs. We focus on comparing the mean optimal vaccination rate in the heterogeneous population with the optimal vaccination rate in the corresponding homogeneous population. We show that there is a critical size for the group with lower relative costs and the mean optimal vaccination in the heterogeneous population is more than in the homogeneous population if and only if the group is larger than the critical size.
Improved hygiene conditions by economic growth and the introduction of the national immunization program for infants have led to variations in hepatitis A antibody prevalence across age groups in Korea. Specifically, individuals in their 20s to 40s have the lowest antibody prevalence. Given that the fatality rate of hepatitis A increases with age, the low immunity level among young adults suggests that, without additional preventive interventions, there is a risk of increased deaths in older age groups in the future. We developed an age-structured transmission model that accounts for age-specific antibody prevalence and fatality rates to assess the impact of adult vaccination, assuming it starts in 2025. We compared vaccination strategies targeting individuals in their 20s to 30s and those in their 40s to 50s, considering that antibody testing costs are incurred for the latter group in Korea. Our study shows that when total costs for vaccination are fixed, vaccinating individuals in their 40s to 50s covers 0.2 times fewer individuals than vaccinating those in their 20s to 30s but reduces deaths by 1.3 to 1.5 times more. When the total vaccine supply is fixed, the total and annual costs of vaccinating individuals in their 40s to 50s are 1.2 times higher than those for the 20s to 30s group, while the reduction in deaths is 1.7 to 1.8 times greater. From the perspective of reducing deaths, vaccinating individuals in their 40s to 50s is more effective than vaccinating those in their 20s to 30s. Furthermore, our research suggests that if an additional vaccination intervention is introduced for individuals in their 20s to 30s, military personnel may continue to receive only a single-dose vaccination, as is currently practiced.
This paper studies a Susceptible-Protected-Infected-Recovered (SPIR) epidemic model incorporating both local and nonlocal diffusion to capture how fear of infection impacts population behaviour. We aim to establish the model’s well-posedness by proving the existence, uniqueness, and positivity of solutions. A key focus is deriving a variational expression for the basic reproduction number mathfrak{R_0}, which acts as a threshold indicator. When mathfrak{R_0} < 1, the disease-free equilibrium is globally stable, meaning the infection dies out. When mathfrak{R_0} > 1, the model predicts disease persistence and the existence of a stable endemic equilibrium, which we prove using Lyapunov functions. We consider two distinct scenarios: one where the diffusion of the susceptible population is absent, and another where the infected population does not diffuse. The model is also compared with the classical SIR model to evaluate the impact of protective measures aimed at reducing mathfrak{R_0} below unity for effective disease control.
Detection and surveillance of emergency animal disease outbreaks play a crucial role in their management. However, monitoring disease spread often relies on farmers self-reporting animal infection. If there are disincentives for farmers to report disease, this delayed surveillance can impact outbreak management outcomes. To understand how farmer reporting rates impact outbreak management strategies, we model the spread of animal disease using an agent-based model. We investigate how varying the reporting rate of infected properties can impact the absolute outcomes of management strategies and which strategy is optimal. Our model considers disease transmission occurring within a property, through animal-to-animal transmission, and between properties, through wind dispersal of fomites and random movement of animals. Using this model, we compare the number of animals culled under four strategies: culling infected properties and animal movement restrictions combined with (1) ring culling, (2) ring testing, (3) ring vaccination (with a perfect vaccine), and (4) ring vaccination (with an imperfect vaccine). Our modelling demonstrates how human behaviour, such as reporting rates, can impact the outcomes from managing emergency animal disease outbreaks. While exact behaviour cannot be predicted for future outbreaks, we can prepare for the next outbreak of emergency animal disease by designing management strategies that are robust to a variety of human behaviours.
We provide an in silico study of stochastic extinction of virus infection. Our work considers a nonspecific antiviral drug that increases the virus clearance rate, and we investigate the effect of this drug on early infection extinction. Synthetic data is generated by a hybrid multiscale framework that applies both continuous and discrete mathematical approaches -- virus spread is described by a partial differential equation, while the cell population is grasped by an agent-based model. The central result we present is the observation, analysis, and explanation of a linear relationship between the virus clearance rate and the probability of early infection extinction. The derivation behind this simple relationship is based on the theory of branching processes.
Stochasticity is introduced to bring new insight into the modelling of population dynamics of diseases. Many systems, in nature, are subject to stochastic perturbations. In this talk, we present differential equations with stochastic perturbations and the updated data estimation method for estimating the transmission rate changes over time. The models for the population dynamics under SEIR epidemic models with stochastic perturbations are presented for the dynamics of the COVID-19 pandemic in Bogota, Colombia. We performed computational experiments to interpret COVID-19 dynamics using actual data for the proposed models. We estimated the model parameters and updated their reported infected and recovered data estimates. (joint work with Andres Rios-Gutierrez )
The Omicron wave of the coronavirus disease 2019 (COVID-19) pandemic posed significant challenges for public health interventions and surveillance efforts due to limited diagnostic capacity to detect all infections. Newfoundland and Labrador (NL), a Canadian province, transitioned between its tiered Alert Level System (ALS) to guide the intensity of non-pharmaceutical interventions (NPIs) as cases surged. However, the true burden of infections and impact of the provincial ALS on Omicron transmission remains unclear due to widespread underreporting driven by restricted access to reverse transcription polymerase chain reaction (RT-PCR) testing. This study resolves the challenge of quantifying Omicron's transmission under conditions of limited testing and surveillance gaps in NL by estimating the true burden of infection using a calibrated mechanistic compartmental model, fit to infection-induced seroprevalence data. By integrating time-varying testing eligibility fractions and estimating transmission parameters, we quantify how the provincial ALS influenced viral transmission between December 15, 2021, and May 26, 2022. Our findings reveal that alert level 4 (ALS-4; the most restrictive under Omicron) was associated with an 85% reduction in Omicron's transmission and marked the only phase in which Omicron incidence decreased while the less restrictive levels rather slowed transmission. The estimated burden of infections (182,534) was over four times the number of reported RT-PCR confirmed cases (41,619), highlighting about 77% substantial under-ascertainment rate. These results contribute to rare empirical evidence that stringent public health restrictions can suppress Omicron transmission in highly vaccinated populations, an outcome that has been observed in only a few jurisdictions and underscores the importance of seroprevalence-informed modelling for policy evaluation.
HIV-1 has nine main subtypes that persist in the infected populations. However, the overall diversity of HIV-1 is much larger thanks to recombination amongst these original subtypes. Recombination has led to over a hundred circulating recombinant forms (CRFS), which may become more prevalent than the two parent strains. We analysed sequences taken from the Los Alamos National Lab HIV database and found different trends for various CRFS depending on the country of origin. However, what drives these differences, is it within-host or between-host dynamics? To investigate, we developed a mathematical model of HIV-1 viral recombination that incorporated within-host dynamics (viral competition and recombination) and between-host dynamics (transmission rates and emergence time in the population) to understand which part of the dynamics matters more. We apply this method to three scenarios regarding CRF emergence in Brazil and China, showing that our model can capture the three scenarios. Given the model assumptions, we also show that withinhost dynamics are an early driver in the emergence of CRFs. However, between-host events will determine the level at which a CRF may be expressed in the population.
Cholera continues to pose a serious risk in developing areas, necessitating strong forecasting to guide public health measures. This research merges qualitative dynamics with machine learning to estimate cholera transmission patterns in Malawi. An epidemic model based on mechanistic principles captures the spread of the disease through parametric calibration. Sensitivity analysis using partial rank correlation coefficients pinpoints critical parameters that affect transmission. The basic reproduction number defines the long-term trends, while bifurcation analysis illustrates how disinfection influences the stability of the disease. To improve predictive accuracy, we combine the mechanistic model with ARIMA and autoregressive neural networks, creating hybrid machine learning models informed by the epidemic context. We generate short-term predictions of cholera cases, showing the advantages of integrating temporal disease dynamics into data-driven approaches. This combined method provides a replicable framework for making forecasts and aids in timely decision-making for epidemic management.
During the COVID-19 pandemic, we saw a widespread adoption of wastewater-based surveillance as a passive, non-invasive tool for monitoring community-level transmission. Early efforts focused on correlating wastewater viral loads with clinical indicators such as reported cases and hospitalizations, often identifying a time lag between wastewater signals and clinical outcomes. As the field evolved, researchers began applying statistical and machine learning models to leverage these signals for short-term forecasting of disease trends. In this study, we contribute to the growing body of work by developing predictive models for SARS-CoV-2, Influenza A, and Respiratory Syncytial Virus (RSV) in British Columbia (BC). Using wastewater viral load as a predictor, we first constructed models to forecast reported case counts. Building on this foundation, we are developing a spatiotemporal hierarchical Bayesian model that integrates reported case data and wastewater signals to predict respiratory-related emergency department (ED) visits across BC. Our work highlights the value of wastewater-based data for early detection and response planning in public health systems.
During the COVID-19 pandemic early years, infection-prevention measures suppressed transmission of seasonal influenza and other respiratory viruses. The early onset and moderate severity of the US 2022-2023 influenza season may have resulted from reduced use of nonpharmaceutical interventions or lower population immunity after two years of limited influenza virus circulation. We used a mathematical model of influenza virus transmission that incorporates vaccine-derived protection against both infection and severe disease, observed hospitalization burden, to estimate the impact of influenza vaccines on healthcare burden. Despite limited data on vaccine effectiveness against infection, our analyses suggest substantial indirect protection, particularly from young adults to other age groups. This is supported by a significant negative correlation between young adult (aged 18-49 years) vaccination rates and observed hospital burden across US states. Assuming reported levels of past vaccine effectiveness against infection and hospitalization, we estimate that influenza vaccines prevented 67,931 [95% confidence interval (CI): 34,182, 95,842] influenza-related hospitalizations nationwide during the 2022-2023 season, with 61% attributable to reduced susceptibility and onward transmission. Among those aged >=65 years, nearly half of averted hospitalizations resulted from vaccinating younger age groups. These findings highlight the need for better estimates of influenza vaccine effectiveness against infection and the potential benefits of increasing young adult influenza vaccination rates to reduce both direct and indirect disease burden.
We study the generation-time distribution, i.e. timing of infection events, in a discrete-time epidemic model with asymptomatic carriers. The progression of the disease is categorized into four phases: the non-infectious latent phase, the infectious asymptomatic phase (a key feature of the model where individuals exhibit mild or no symptoms), the infectious symptomatic phase, and lastly, the immune phase. We introduce a versatile non-Markovian system with generic waiting times at infected stages and transmission rates depending on the elapsed time since infection. The basic reproduction number is derived from a renewal equation for the (sequence of) asymptomatic hosts, whose expression gives the probability distribution of the time between new cases in a chain of infection transmission (generation time). For illustration purposes, we consider Weibull distributions which include both geometric and fixed-length distributions as particular cases. For memoryless waiting times (geometric distribution), we have investigated the evolution of infection transmission before and after symptom onset. Given that individuals can develop symptoms and die from the disease, we consider disease-induced mortality as a measure of virulence and assume it is positively correlated with a weighted average transmission rate. Our findings indicate that the infection transmission rate is consistently higher during the symptomatic phase. However, in some scenarios, the majority of infections occur prior to symptom onset. [1] A Discrete Model for the Evolution of Infection Prior to Symptom Onset, J. Ripoll and J. Font, Mathematics 2023, 11(5), 1092. [2] Generation time in a discrete epidemic model with asymptomatic carriers: beyond geometric waiting-times, J. Ripoll and J. Saldaña, IMAE-preprint.
Self-isolation and stay-at-home measures are crucial for curbing the spread of contagious pathogens while vaccines are being developed. Furthermore, research during the 2019-22 coronavirus pandemic (COVID-19) emphasizes that proper enforcement and timely lifting of these measures are vital for effective disease management. In this context, we analyzed a simple dynamical system to understand how an epidemic progresses by isolating susceptible individuals (confinement) and reintroducing them to infection (deconfinement). This model captures the overall magnitude and direction of flows between confined and deconfined groups—akin to osmosis—leading to a dimensionless quantity defined as confinement tonicity. Our mathematical analysis suggests that confinement tonicity influences the final epidemic size, providing insights into careful quarantine management for effective disease control.
The accrual of nucleotide substitutions in pathogen genomes accompanies their transmission through host populations. Because lineages with higher fitness tend to transmit rapidly to new hosts before incurring very many substitutions, large numbers of related sequences are usually interpreted as evidence of transmission success. Quantities like the local branching index (LBI) aim to identify successful lineages in this way by scoring sequences according to the number of close relatives captured in the dataset. While statistics like LBI are easily calculated from a given phylogenetic tree (or a distribution of trees), observation errors related to sampling bias and censoring may introduce spurious signals of transmission success. To disentangle these effects, we use stochastic compartmental models to simulate outbreaks and generate distributions of phylogenies under a variety of testing programs (such as surveillance of symptomatic cases, or cross-sectional prevalence studies). By characterizing the types of phylogenies expected under these situations, we can work towards a clearer understanding of the types of signals that are likely to be detected with sequence data.
Dengue fever has emerged as an increasingly alarming public health challenge, further complicated by the impacts of climate change on control efforts. Yet, the full extent of climate's impact on dengue incidence remains poorly understood. To investigate this, we employed an advanced causal inference method to 16 regions in the Philippines, selected for their diverse climatic conditions. Unlike previous methods for detecting regulatory relationships, this method is capable of detecting nonlinear and joint effects of temperature and rainfall to dengue incidence. We found that temperature consistently increased dengue incidence throughout all the regions, while rainfall effects differed depending on location. Further analysis showed that this pattern is due to the variation in dry season length, a factor previously overlooked. Specifically, our results showed that regions with low variation in dry season length experience a negative impact of rainfall on dengue incidence likely due to strong flushing effect on mosquito habitats, while regions with high variation in dry season length experience a positive impact, likely due to increased mosquito breeding sites. This study offers a fresh perspective on the relationship between climate and dengue incidence, emphasizing the need for tailored prevention strategies based on local climate conditions.
This study investigates the spatial-temporal heterogeneity in the relationship between human mobility and COVID-19 transmission across 229 regions in South Korea during six epidemic waves from January 2020 to September 2022. While previous research primarily focused on the early stages of the pandemic and the impacts of mobility restrictions, our study utilizes mobility data from SK Telecom and COVID-19 case data from the Korea Disease Control and Prevention Agency to provide a more comprehensive analysis. We applied empirical mode decomposition (EMD) and clustering analysis to classify regional mobility patterns and conducted cross-correlation analysis to assess the relationship between mobility and confirmed cases. The findings indicate that incoming mobility significantly influenced the number of confirmed cases in urban and densely populated areas, whereas rural regions exhibited contrasting patterns. Moreover, these relationships evolved across different epidemic waves, highlighting the influence of regional characteristics and public health interventions. This study underscores the need to consider spatial-temporal heterogeneity in mobility-transmission dynamics to develop tailored public health strategies and enhance preparedness for future pandemics.
Forecasting the occurrence and absence of novel disease outbreaks is essential for disease management, yet existing methods are often context-specific, require a long preparation time, and non-outbreak prediction remains understudied. To address this gap, we propose a novel framework using a feature-based time series classification (TSC) method to forecast outbreaks and non-outbreaks. We tested our methods on synthetic data from a Susceptible–Infected–Recovered (SIR) model for slowly changing, noisy disease dynamics. Outbreak sequences give a transcritical bifurcation within a specified future time window, whereas non-outbreak (null bifurcation) sequences do not. We identified incipient differences, reflected in 22 statistical features and 5 early warning signal indicators, in time series of infectives leading to future outbreaks and non-outbreaks. Classifier performance, given by the area under the receiver-operating curve (AUC), ranged from 0.99 for large expanding windows of training data to 0.7 for small rolling windows. The framework is further evaluated on four empirical datasets: COVID-19 incidence data from Singapore, 18 other countries, and Edmonton, Canada, as well as SARS data from Hong Kong, with two classifiers exhibiting consistently high accuracy. Our results highlight detectable statistical features distinguishing outbreak and non-outbreak sequences well before potential occurrence, in both synthetic and real-world datasets presented in this study.
The increase in Plasmodium vivax malaria cases in Korea highlights the necessity to reevaluate intervention strategies as climate patterns change. In 2024, confirmed cases rose by 37% compared to the previous three years' average, along with an increase in vector mosquito populations. In response, the Korea Disease Control and Prevention Agency (KCDA) expanded designated malaria risk areas in Seoul. Effective control depends on optimizing non-pharmaceutical interventions with primaquine-based treatment. As tafenoquine emerges as a potential alternative treatment, evaluating its impact on malaria transmission, relapse rate and cost-effectiveness within public health systems is essential. To address these issues, we developed a mathematical model incorporating climate variability to assess the effectiveness of non-pharmaceutical interventions under different climate scenarios. Using the Improved Multi-Objective Differential Evolution (IMODE) algorithm, we analyzed the optimal interventions based on observed malaria control measures. Our results suggest that optimal intervention strategies can significantly reduce malaria transmission and relapse rate, highlighting the cost-effectiveness of tafenoquine and optimal intervention approaches in Korea’s malaria control measures.
Abstract Background: Healthcare accessibility is a key determinant of health outcomes during pandemics. Disparities in access may contribute to indirect excess mortality beyond reported COVID-19 deaths. This study quantitatively examines the impact of regional healthcare accessibility on non-COVID excess mortality in South Korea using a mathematical modeling approach. Methods: We first performed hierarchical clustering based on the average travel time to various healthcare facilities, classifying regions into two groups: Cluster 0 (high accessibility) and Cluster 1 (low accessibility). A CatBoost model trained on 2014–2019 data predicted expected deaths for 2020–2022, and excess mortality was calculated as the difference between observed and predicted deaths. Finally, multiple linear regression was then used to evaluate the association between accessibility time and non-COVID excess mortality. Results: Our analysis revealed that regions with high healthcare accessibility (Cluster 0) exhibited excess mortality patterns that closely aligned with reported COVID-19 deaths. In contrast, regions with lower accessibility (Cluster 1) experienced a significant increase in non-COVID excess mortality, particularly during the Omicron surge (fifth and sixth pandemic waves). The regression analysis demonstrated that longer healthcare accessibility times were significantly associated with higher non-COVID excess mortality in later pandemic stages. Conclusion: This study demonstrates that regional disparities in healthcare accessibility contribute to indirect excess mortality during pandemics. The findings highlight the importance of targeted policy interventions, such as strengthening healthcare infrastructure and expanding telemedicine, to reduce health inequalities and enhance public health resilience in future crises.
Repeating patterns, such as hair follicles and bristles play important roles in the lives of animals. These structures help animals to optimally sense their environment. Notch signaling is known to control these patterns. Primarily, Notch signaling is a local communication between neighboring cells in contact (signal-sending and signal-receiving cells). The local communication between cells in contact is not able to explain all the complex biological patterns observed. Further studies reveal long-range communication between cells using actin-based filopodia called cytonemes. The precise understanding of how the dynamics of filopodia such as protrusion and retraction lead to notch-delta activation remains unclear. In this work, we develop a mathematical model to help unravel the mystery of this long-range communication between cells. Student: William Ebo Annan Advisors: Prof. Emmanuel O.A. Asante & Prof. Ginger Hunter.
Gene regulation, affected by random molecular fluctuations, is often modeled assuming DNA is evenly distributed in the nucleus—an unrealistic simplification. We found that when key molecules move slowly, these models fail unless uneven spatial distribution is included, which slows simulations. We explored simplification techniques to speed up the process while keeping accuracy. This study stresses the need for tools balancing efficiency and precision in modeling gene regulation with spatial differences.
Systems of intracellular biochemical reactions are complex, often involving components that cannot be directly measured. Representing these systems as networks, with nodes representing biochemical species and edges their reactions, helps quantitatively characterize their function and effects of dysregulation. Causal discovery methods can uncover functional interactions within these networks from purely observational data, detecting hidden effects from partial observations. These effects appear as common causes of observed variables, or through time-lagged effects from intermediate causes. We benchmark the causal discovery method temporal Multivariate Information-based Inductive Causation (tMIIC) alongside other state-of-the-art tools, for time series data from biochemical kinetic models. Our results demonstrate tMIIC’s high recall in identifying interactions within toy reaction networks. By selectively omitting data, we consider both latent confounders (the standard choice for benchmarking these methods) and unobserved species participating in reactions. tMIIC detects latent confounders using bidirected edges, and unobserved species through time-delayed edges, locating hidden effects and estimating their typical timescales. Finally, we extend these benchmarks to reconstruct an experimentally calibrated model of the epidermal growth factor receptor signalling network – a system frequently dysregulated in cancer. Altogether, our work showcases the feasibility and usefulness of causal discovery methods like tMIIC for data-driven mathematical modelling of biochemical reactions.
Sickle cell disease (SCD) is a group of inherited health conditions that affect the red blood cells. Millions across the globe are affected by SCD. More than 100,000 Americans and nine out of ten people in the United States who have SCD are of African descent. Individuals that carry SCD produce these abnormally shaped red blood cells (RBC) which can adversely affect the body. These red blood cells shaped like sickles do not live as long as healthy RBCs and can cause blockages in blood vessels that can lead to pain. Managing pain episodes is the focus of our current research. Valrie et al. (2021) showed that there were correlations between sleep quality and pain the next day and also positive affect and pain. Positive affect (PA) is measured using self-reported scales to evaluate the level of positive emotions a person is feeling at a certain time. Therefore, we have developed different mathematical models to study sickle cell disease pain, one of which shows the relationship between sleep and pain and the other that focuses on positive affect. Previous studies have considered PA as a mediator between sleep and pain, however for our research, we will treat positive affect as a potential driver to predict pain. This work utilizes diary data from pediatric patients and considers differences in the adolescent and children subpopulations.
Interacting particle systems, also known as agent-based models (ABMs), represent one category of dynamical systems that are used to study a wide range of physical phenomena across multiple scales. Examples from science and engineering include cell migration, swarm robotics, social psychology, and animal migration patterns and interactions. A ubiquitous feature of such systems is that they exhibit a form of emergence: local interactions leading to large-scale coordination. A fundamental scientific question is thus to understand the local interactions that give rise to the observed emergent dynamics. We are interested in methods for learning interactions generally, which can describe any ABM defined by an interaction kernel without making any additional assumptions about the analytical form of this kernel (i.e. it is non-parametric). The advantage of this kernel-based approach is that it incorporates the underlying physics of the model (i.e. collective dynamics), which more general approaches may ignore, potentially limiting their effectiveness. We propose to extend a non-parametric statistical learning approach for learning the interaction kernel for systems with both self-propulsion and collective dynamics, given an observed set of trajectories. First, we parametrically learn the intra-agent force while simultaneously inferring the interaction kernel non-parametrically. The method is validated on two well-known models. We extended this approach to learn the intra-agent force non-parametrically. Also, we explored how to identify the best-fit model among all possible variations for learning interacting particle collective motion based on observations.. Also, we will introduce an alternative neural network framework to the existing non-parametric statistical learning approach.
As wildfires increase in frequency and intensity due to climate change, so does the need to create better forecasts. The wildfire perimeter is seldom fully observable, but the geospatial locations of first responder and 911 civilian cellphone calls provide a sparse representation of the wildfire in real time. We use these calls to estimate the complete fire perimeter at specific times. We present a state estimation method that yields smaller reconstruction errors than existing methods. Using the reconstruction as an initial condition, we run a cellular automaton to predict the future state of the fire. As a case study, we use calls from Maui during the devastating wildfires in August 2023 to predict the final fire perimeter.
One of the difficulties in eliminating Plasmodium vivax malaria lies in its ability to cause recurrent infections following the activation of dormant parasites (relapse). However, this can be confused with recurrent infections due to treatment failure (recrudescence), or a new infectious mosquito bite (reinfection). Distinguishing the cause of recurrent Plasmodium vivax malaria in each patient is critical for malaria control efforts, such as efficacy studies of drug treatments. We address this need by developing a statistical tool to infer the cause of Plasmodium vivax recurrent malaria from genetic data, implemented through the R package Pv3Rs. Each mode of malaria recurrence – relapse, recrudescence, and relapse – feature different levels of genetic relatedness between parasites. We use Bayesian hierarchical modelling to translate genetic relatedness in observed data to interpretable probabilities for each mode of malaria recurrence. We illustrate the utility of our model by applying it to Plasmodium vivax microsatellite marker data of acute malaria patients treated with high-dose primaquine. The ability to probabilistically resolve the cause of recurrent malaria helps provide more accurate failure rates of drug treatments.
Genetically engineered bacteria are increasingly utilized to manufacture products that are difficult, expensive, or impractical to synthesize chemically. These products have potential applications ranging from medicine to sustainability. However, metabolic pathway introduction, extensive feedback mechanisms in the cell, and evolutionary forces complicate the engineering of bacterial strains that are well-suited for the task. We need tools that will enable us to anticipate these challenges, as well as increase efficiency and enable novel design. A recently published large-scale model of Escherichia coli has enabled us to simulate many distinct cellular processes and capture their complex interactions on a system-wide level. This model incorporates decades of heterogeneous data collection from E. coli literature to fit over 19,000 parameters for the mechanistic ordinary differential equations describing processes in the cell. We now introduce components related to genetic engineering, with an initial focus on chromosome modification. In this poster, we describe preliminary work analyzing the trade-offs between maximizing exogenous protein production and preserving cell health. Our numerical experiments varying the expression level of a single gfp gene reveal how exogenous gene products sequester resources in key cellular processes. We anticipate that these methods will set the stage for large-scale computational genetic engineering design tools as they develop and expand.
The Svalbard reindeer is a subspecies of Rangifer tarandus that is endemic to the arctic island, making them vulnerable to the impacts of climate change. Among the semi-isolated coastal populations, juvenile dispersal is crucial to maintaining viability, particularly with more rain-on-snow events making it challenging to access food. However, the absence of sea ice is a potential barrier to dispersal. Using step selection functions, we can understand how parental mimicry and individual exploration conditioned on habitat covariates drive dispersal among juveniles. Incorporating learning into step selection functions is an emerging area of research, allowing for a deeper understanding of animal movement behaviour. This project will develop new techniques for step selection functions, as well as providing key insights regarding learning and dispersal of juvenile Svalbard reindeer.
Repeating patterns, such as hair follicles and bristles play important roles in the lives of animals. These structures help animals to optimally sense their environment. Notch signaling is known to control these patterns. Primarily, Notch signaling is a local communication between neighboring cells in contact (signal-sending and signal-receiving cells). The local communication between cells in contact is not able to explain all the complex biological patterns observed. Further studies reveal long-range communication between cells using actin-based filopodia called cytonemes. The precise understanding of how the dynamics of filopodia such as protrusion and retraction lead to notch-delta activation remains unclear. In this work, we develop a mathematical model to help unravel the mystery of this long-range communication between cells. Student: William Ebo Annan Advisors: Prof. Emmanuel O.A. Asante & Prof. Ginger Hunter.
Modeling with systems of differential equations requires prior knowledge to create a fully specified model reflecting our understanding of biological systems, but sometimes we don’t have a complete understanding of the systems we are interested in. If we have time series data of our variables of interest, multilayer perceptron models can take the place of unknown terms in our equations to produce solutions that fit the data well but due to the nature of multilayer perceptron models are not very interpretable. Interpretability can be improved by combining the simpler compositionality of Kolmogorov Arnold Networks with symbolic regression, allowing for the discovery of unknown terms from time series data. In this poster, we leverage the interpretability of Kolmogorov Arnold Networks with symbolic regression to recover a logistic growth term from time series data generated by a predator-prey model.
Accurate ecological forecasts provide useful insights to inform policy and management, but building models to produce these forecasts is challenging. Modelling approaches can vary from mechanistic models that attempt to capture the underlying ecological processes to purely phenomenological or statistical models that rely on inferences from data. These different approaches are likely to have different strengths depending on the metric being predicted, the amount of data for training, and the time horizon of the prediction. In particular, too strong reliance on past data may lead to incorrect inferences about the future of ecosystems under novel conditions, such as those induced by climate change and anthropogenic disturbances. We here study several models of different paradigms, including neutral models, to predict Mountain pine beetle (MPB) infestations in Alberta, Canada. MPB life history makes mathematical modeling challenging, as they use complicated chemical signalling processes and occasionally disperse over very long distances above the tree canopy. During a recent hyperepidemic in neighbouring British Columbia, MPB were able to overcome the natural border of the Rocky Mountains and spread into Alberta. Alberta dedicated extensive resources to monitor and control this spread, including helicopter surveys of infested trees. We use this data to study the predictive accuracy of several models that range in complexity and mechanistic basis. We discuss general trends in model performance with the aim of providing practice advice about the types of models that may achieve the greatest predictive accuracy given different data availability, target year, and forecast horizons.
Invasive species pose a significant threat to biodiversity, the environment, and the economy. They are expensive to manage and monitor, in the UK alone the estimated annual cost to the economy is £4 billion. The spread of invasive species is increasing at unprecedented rates, as a result of expanding human trade networks and climate change. One invasive species of note is the oak processionary moth (OPM), which became established in the UK in 2006 through accidental importation. OPMs are harmful defoliators of oak trees, leaving them vulnerable to other stressors and diseases. They are also harmful to humans; the caterpillars have urticating hairs which can cause breathing difficulties. The eradication of OPM in the UK has been deemed unfeasible with the current management strategy focused on containing their spread. In partnership with Fera Science, we are combining mathematical and statistical models to describe and predict the spread of OPM across the UK and to inform future management strategies. This poster will showcase our work on an agent-based (individual-based) modelling approach for capturing OPM spread across the South-East of England. This model uses a lattice-based grid where a cell is either infested or susceptible (analogous to an SI model) to OPM, with cells becoming infested based on their distance to infested cells under the assumption of different dispersal kernels. We can then use the model to guide new management strategies and scenario test which may allow better containment.
Immune cells observed under a microscope often exhibit motion that deviates from simple random (Brownian) walks, yet the precise factors driving these deviations remain poorly understood. Statistical approaches, such as hidden Markov models, segment cell trajectories into multiple pure Brownian motion regimes and estimate a diffusion coefficient for each. While these methods provide insight into short-term movement changes and can be used to predict future positions, they offer limited explanation of the underlying biological motivation. Here, we propose a novel framework based on a transport equation to probe the fundamental causes of anomalous diffusion in immune cell trajectories. By fitting a generalized model to data, we can distinguish whether deviations from pure diffusion arise from a correlation in the cell’s motion—implying a “memory” of previous steps—or from a bias driven by an external factor, such as a cytokine gradient. In the latter scenario, immune cells may adjust their paths when approaching target cells (e.g., cancer cells), giving rise to a biased random walk. Moreover, the goal is to come up with a framework that accomodates the possibility of both memory effects and bias (a biased correlated random walk), enabling a more comprehensive description of immune cell motility. From our initial tests, it looks like this transport-based approach can spot unique signs of correlation or bias in immune cell movement just by analyzing time-lapse data. This could help us better understand how immune cells navigate their environments and may even open new avenues for guiding their behavior in therapeutic settings.
Enzyme inhibition analysis is essential in drug development and food processing, necessitating precise estimation of inhibition constants. Traditionally, these constants are estimated through experiments using multiple substrate and inhibitor concentrations, but inconsistencies across studies highlight a need for a more systematic approach to set experimental designs across all types of enzyme inhibition. Here, we addressed this by analyzing the error landscape of estimations in various experimental designs. We found that nearly half of the conventional data is dispensable and even introduces bias. Instead, by incorporating the relationship between IC50 and inhibition constants into the fitting process, we found that using a single inhibitor concentration greater than IC50 suffices for precise estimation. This novel IC50-based optimal approach, which we name 50-BOA, substantially reduces (>75%) the number of experiments required while ensuring precision and accuracy. Additionally, we provide a user-friendly package that implements the 50-BOA.
Biological condensates are membraneless organelles within the cell or the nucleus which perform an array of different tasks and typically consist of DNA/RNA and protein. Actin is a protein that exists in most eukaryotic cells and transitions between monomeric and filamentous states. In its filamentous state, actin forms networks, which perform vital tasks inside the cell. Recent research has shown interactions between biological condensates and cytoskeletal filaments, such as actin. The focus of these works was on morphological changes of condensates, transportation of condensates along cytoskeletal structures and on bundling of actin filaments inside condensates. While works have shown that condensates facilitate actin polymerization, a theoretical mathematical description of the cooperativity between condensates and actin polymerization is still missing. In this work we use a master equation to capture the polymerization kinetics of actin in a multicompartment system of condensates and dilute phase containing monomeric and filamentous actin. We believe that the model will allow us to make predictions about the number and length of fibers polymerized inside condensates. The findings of this study will hopefully further our understanding of the cooperative behavior between actin and biological condensates and help us understand how biological condensates are involved in the creation and maintenance of the actin networks inside the cell.
Inspired by the dynamics of Influenza A attachment to the epithelial cells of the upper respiratory tract, we are developing a dynamic biophysical model of virus attachment to cell surface receptors. Major challenges to modelling include modelling multiple ligands distributed across the viral surface, and developing a dynamic binding and unbinding model which incorporates forces applied to the ligand-receptor pair. In this presentation, I will describe some of these challenges and explain how our computational approach based on the principles of Steven Andrews SMOLDYN is being used to overcome them. I will also describe how we envisage using our framework to study practical questions about how viruses penetrate mucosal and ciliary layers in the upper respiratory tract, and how this framework can be extended to virus-like systems such as engineered nanoparticles for drug delivery.
Safta et al (J.~Comp.~Phys. 2015) introduced a hybrid discrete/continuous representation of Kolmogorov's forward operator, $mathcal{L}$, for numerically simulating the evolution of probability distributions on state spaces spanning both large and small numbers of molecules. Motivated by first-passage-time (e.g. extinction time) and exit-distribution problems, we extend Safta et al's approach to establish a hybrid discrete/continuum representation of Kolmogorov's backward operator $mathcal{L}^dagger$, the formal adjoint of $mathcal{L}$ also known as the Markov process generator, or the stochastic Koopman operator. We apply our coarsened backward operators to several birth-death processes of increasing complexity, leveraging exact results where available to evaluate their speed and accuracy.
Previous studies have demonstrated that sleep patterns, when combined with mathematical models describing sleep homeostasis and circadian rhythms, can predict individual alertness and mood. However, these studies have primarily focused on adults, leaving it unclear whether these approaches are applicable to adolescents. In this study, we investigate whether adolescents' alertness and mood can be predicted solely from their sleep patterns.
Memory plays a crucial role in shaping animal movement behaviors. There are powerful new models emerging for understanding memory effects via continuous-time stochastic differential equations (SDEs). However, these models are non-Markovian and computationally difficult to implement as a result of intractable likelihood functions. To address such issues, we establish a connection between discrete-time and continuous-time animal movement models and introduce a conditional Markov process, facilitated by a so-called filtration or dynamically updated cognitive map. Such cognitive maps are allocentric mental representations of an individual’s surroundings, which change over time, and have been shown to influence animal movement and habitat use. We reformulate the non-Markovian memory-mediated animal movement SDE (M3) model originally developed by Fagan et al. (2023) by incorporating cognitive maps to create a Markov process conditioned on the associated cognitive map. This new formulation facilitates parameter inference and enhances computational efficiency. Introducing cognitive maps into the M3 model transforms the non-Markovian SDE model to a conditionally Markovian model that is coupled to an auxiliary filtration that is described as a cognitive map. This cognitive map, in turn, is described by a spatially distributed system of ordinary differential equations, which describe how the animal's spatial memory changes over time. The reformulation offers practical advantages by making the likelihood function computationally tractable and avoiding large-scale integrals. Simulations show that our cognitive map-based model (CM-M3) preserves the key movement patterns of the original M3 model—bounded wandering, convergence, and cyclic paths. We also demonstrate improved simulation speeds using optimized cell sizes and interpolation, and explore parameter inference methods that address the issue of intractable likelihoods.
Gene regulation, affected by random molecular fluctuations, is often modeled assuming DNA is evenly distributed in the nucleus—an unrealistic simplification. We found that when key molecules move slowly, these models fail unless uneven spatial distribution is included, which slows simulations. We explored simplification techniques to speed up the process while keeping accuracy. This study stresses the need for tools balancing efficiency and precision in modeling gene regulation with spatial differences.
Organoids are powerful models for studying cellular self-organisation and tissue morphogenesis, with applications in disease modelling, personalized medicine, and drug screening. However, their growth and development mechanisms remain incompletely understood, contributing to variability in lab-grown organoids, and challenging the design of new microdevices to reliably represent the complexity of in-vivo conditions. This poster presents recent advances in the development of mathematical modelling tools to better understand organoid behaviour and support the development of organoid-on-a-chip systems. We introduce two new agent-based models developed using the multi-scale simulation framework Chaste (https://chaste.github.io/): (1) a 2D organoid model incorporating updated cell-cycle dynamics based on in vitro observations, and (2) a 3D organoid-on-a-chip model that integrates the organoid with its microfluidic environment. Both models are evaluated for their ability to reproduce budding structures—a key morphological feature of intestinal organoids. To support model validation, we developed a novel machine learning algorithm that automates the counting of budding structures in both experimental and simulated images. Our results demonstrate that the in silico models can replicate budding counts observed in our new in vitro data. These tools enhance the identification and quantification of key morphological features, enabling deeper comparisons between computational models and laboratory data, with potential benefits for experimental and mathematical model design and interpretation.
Systems of intracellular biochemical reactions are complex, often involving components that cannot be directly measured. Representing these systems as networks, with nodes representing biochemical species and edges their reactions, helps quantitatively characterize their function and effects of dysregulation. Causal discovery methods can uncover functional interactions within these networks from purely observational data, detecting hidden effects from partial observations. These effects appear as common causes of observed variables, or through time-lagged effects from intermediate causes. We benchmark the causal discovery method temporal Multivariate Information-based Inductive Causation (tMIIC) alongside other state-of-the-art tools, for time series data from biochemical kinetic models. Our results demonstrate tMIIC’s high recall in identifying interactions within toy reaction networks. By selectively omitting data, we consider both latent confounders (the standard choice for benchmarking these methods) and unobserved species participating in reactions. tMIIC detects latent confounders using bidirected edges, and unobserved species through time-delayed edges, locating hidden effects and estimating their typical timescales. Finally, we extend these benchmarks to reconstruct an experimentally calibrated model of the epidermal growth factor receptor signalling network – a system frequently dysregulated in cancer. Altogether, our work showcases the feasibility and usefulness of causal discovery methods like tMIIC for data-driven mathematical modelling of biochemical reactions.
Sickle cell disease (SCD) is a group of inherited health conditions that affect the red blood cells. Millions across the globe are affected by SCD. More than 100,000 Americans and nine out of ten people in the United States who have SCD are of African descent. Individuals that carry SCD produce these abnormally shaped red blood cells (RBC) which can adversely affect the body. These red blood cells shaped like sickles do not live as long as healthy RBCs and can cause blockages in blood vessels that can lead to pain. Managing pain episodes is the focus of our current research. Valrie et al. (2021) showed that there were correlations between sleep quality and pain the next day and also positive affect and pain. Positive affect (PA) is measured using self-reported scales to evaluate the level of positive emotions a person is feeling at a certain time. Therefore, we have developed different mathematical models to study sickle cell disease pain, one of which shows the relationship between sleep and pain and the other that focuses on positive affect. Previous studies have considered PA as a mediator between sleep and pain, however for our research, we will treat positive affect as a potential driver to predict pain. This work utilizes diary data from pediatric patients and considers differences in the adolescent and children subpopulations.
Interacting particle systems, also known as agent-based models (ABMs), represent one category of dynamical systems that are used to study a wide range of physical phenomena across multiple scales. Examples from science and engineering include cell migration, swarm robotics, social psychology, and animal migration patterns and interactions. A ubiquitous feature of such systems is that they exhibit a form of emergence: local interactions leading to large-scale coordination. A fundamental scientific question is thus to understand the local interactions that give rise to the observed emergent dynamics. We are interested in methods for learning interactions generally, which can describe any ABM defined by an interaction kernel without making any additional assumptions about the analytical form of this kernel (i.e. it is non-parametric). The advantage of this kernel-based approach is that it incorporates the underlying physics of the model (i.e. collective dynamics), which more general approaches may ignore, potentially limiting their effectiveness. We propose to extend a non-parametric statistical learning approach for learning the interaction kernel for systems with both self-propulsion and collective dynamics, given an observed set of trajectories. First, we parametrically learn the intra-agent force while simultaneously inferring the interaction kernel non-parametrically. The method is validated on two well-known models. We extended this approach to learn the intra-agent force non-parametrically. Also, we explored how to identify the best-fit model among all possible variations for learning interacting particle collective motion based on observations.. Also, we will introduce an alternative neural network framework to the existing non-parametric statistical learning approach.
As wildfires increase in frequency and intensity due to climate change, so does the need to create better forecasts. The wildfire perimeter is seldom fully observable, but the geospatial locations of first responder and 911 civilian cellphone calls provide a sparse representation of the wildfire in real time. We use these calls to estimate the complete fire perimeter at specific times. We present a state estimation method that yields smaller reconstruction errors than existing methods. Using the reconstruction as an initial condition, we run a cellular automaton to predict the future state of the fire. As a case study, we use calls from Maui during the devastating wildfires in August 2023 to predict the final fire perimeter.
One of the difficulties in eliminating Plasmodium vivax malaria lies in its ability to cause recurrent infections following the activation of dormant parasites (relapse). However, this can be confused with recurrent infections due to treatment failure (recrudescence), or a new infectious mosquito bite (reinfection). Distinguishing the cause of recurrent Plasmodium vivax malaria in each patient is critical for malaria control efforts, such as efficacy studies of drug treatments. We address this need by developing a statistical tool to infer the cause of Plasmodium vivax recurrent malaria from genetic data, implemented through the R package Pv3Rs. Each mode of malaria recurrence – relapse, recrudescence, and relapse – feature different levels of genetic relatedness between parasites. We use Bayesian hierarchical modelling to translate genetic relatedness in observed data to interpretable probabilities for each mode of malaria recurrence. We illustrate the utility of our model by applying it to Plasmodium vivax microsatellite marker data of acute malaria patients treated with high-dose primaquine. The ability to probabilistically resolve the cause of recurrent malaria helps provide more accurate failure rates of drug treatments.
Genetically engineered bacteria are increasingly utilized to manufacture products that are difficult, expensive, or impractical to synthesize chemically. These products have potential applications ranging from medicine to sustainability. However, metabolic pathway introduction, extensive feedback mechanisms in the cell, and evolutionary forces complicate the engineering of bacterial strains that are well-suited for the task. We need tools that will enable us to anticipate these challenges, as well as increase efficiency and enable novel design. A recently published large-scale model of Escherichia coli has enabled us to simulate many distinct cellular processes and capture their complex interactions on a system-wide level. This model incorporates decades of heterogeneous data collection from E. coli literature to fit over 19,000 parameters for the mechanistic ordinary differential equations describing processes in the cell. We now introduce components related to genetic engineering, with an initial focus on chromosome modification. In this poster, we describe preliminary work analyzing the trade-offs between maximizing exogenous protein production and preserving cell health. Our numerical experiments varying the expression level of a single gfp gene reveal how exogenous gene products sequester resources in key cellular processes. We anticipate that these methods will set the stage for large-scale computational genetic engineering design tools as they develop and expand.
The Svalbard reindeer is a subspecies of Rangifer tarandus that is endemic to the arctic island, making them vulnerable to the impacts of climate change. Among the semi-isolated coastal populations, juvenile dispersal is crucial to maintaining viability, particularly with more rain-on-snow events making it challenging to access food. However, the absence of sea ice is a potential barrier to dispersal. Using step selection functions, we can understand how parental mimicry and individual exploration conditioned on habitat covariates drive dispersal among juveniles. Incorporating learning into step selection functions is an emerging area of research, allowing for a deeper understanding of animal movement behaviour. This project will develop new techniques for step selection functions, as well as providing key insights regarding learning and dispersal of juvenile Svalbard reindeer.
This study investigates the finite-time synchronization of complex-valued fractional-order memristive neural networks (CVFOMNNs) with time-varying delays. By separating the model into real and imaginary parts, we design a unified sliding-mode surface and construct a suitable sliding-mode controller to synchronize the drive system state trajectories to those of the response system. We prove that the novel super-twisting sliding mode controller forces the error system to reach a sliding surface in a finite time and decreases the chattering effect. Moreover, the LMI conditions are formulated to guarantee the finite-time synchronization of the CVFOMNNs with time-varying delays by using a less conservative Lyapunov-Krasovskii function. Finally, numerical simulations are presented to validate the effectiveness of the theoretical results.
Childhood absence epilepsy (CAE) is a pediatric generalized epilepsy disorder characterized by brief episodes of impaired consciousness and distinctive 2.5--5 Hz spike-wave discharges (SWDs) on electroencephalography. Although CAE often remits spontaneously during adolescence, the mechanisms driving remission remain poorly understood to effectively guide early intervention practices. Progesterone and its neuroactive metabolite allopregnanolone (ALLO) have been implicated in modulating absence seizure activity. Using a conductance-based thalamocortical model, we previously demonstrated that ALLO enhances GABAa receptor-mediated inhibition, resolving SWDs and supporting the hypothesis that pubertal hormonal shifts may facilitate remission. However, not all patients experience remission despite similar hormonal changes. To investigate mechanisms of resistance to ALLO-mediated remission, we developed an enhanced thalamocortical model that incorporates a layered cortical structure to examine regional cortical heterogeneity and frontocortical connectivity. Our results suggest that non-resolving CAE may result not only from increased frontocortical connectivity but also from the underlying cellular composition of the network. In particular, a higher proportion of bursting-type neurons may prevent the therapeutic effects of allopregnanolone. This work highlights the role of network-level properties in influencing disease outcomes and demonstrates the utility of computational modelling in exploring divergent disease trajectories where empirical models remain limited.
This study investigates the finite-time synchronization of complex-valued fractional-order memristive neural networks (CVFOMNNs) with time-varying delays. By separating the model into real and imaginary parts, we design a unified sliding-mode surface and construct a suitable sliding-mode controller to synchronize the drive system state trajectories to those of the response system. We prove that the novel super-twisting sliding mode controller forces the error system to reach a sliding surface in a finite time and decreases the chattering effect. Moreover, the LMI conditions are formulated to guarantee the finite-time synchronization of the CVFOMNNs with time-varying delays by using a less conservative Lyapunov-Krasovskii function. Finally, numerical simulations are presented to validate the effectiveness of the theoretical results.
Synthetic receptors are engineered proteins with potential applications to studying cell-cell interactions and cancer cell therapy. One promising research direction is engineering the Notch receptor, a transmembrane protein that can detect extracellular signals such as antigens or other ligands, and convert them to intracellular signals to activate expression of certain genes. Both the intracellular and extracellular domains can be engineered and replaced with alternative domains, creating the family of modified Notch receptors known as synthetic Notch (synNotch). SynNotch has the potential to improve chimeric antigen receptor (CAR) T-cell therapy by tuning binding affinity to a specific cancer antigen and minimizing off-target effects. We propose an ordinary differential equation model of synNotch receptor activity that has predictive value of how custom cell response behaviors can be programmed. The mathematical model couples activation dynamics on a fast timescale, characteristic of receptor-ligand interactions, and of slower downstream gene expression dynamics. Local and global sensitivity analyses indicate model parameters that yield the greatest variability in downstream results, indicating their potential to be engineered for different functions. Specifically, we find that ligand association and ligand-dependent activation have the greatest potential for modulating transcription factor release.
Circulating tumor DNA (ctDNA) has emerged as a promising biomarker for monitoring cancer progression and treatment response in real time. In anal squamous cell carcinoma (ASCC), where 80-90% of cases are linked to human papillomavirus (HPV), ctDNA demonstrates high sensitivity in tracking disease dynamics, often detecting progression earlier than imaging while enabling frequent assessment and correlating with tumor burden. Our study examined how patient-specific modeling of ctDNA dynamics can predict time to progression in HPV-associated ASCC. We analyzed longitudinal data from 32 ASCC patients receiving immunotherapy every 3 weeks for up to 2 years, exploring correlations between tumor volume and ctDNA levels. We developed a mathematical model calibrated to patient-specific tumor volume and ctDNA dynamics during immunotherapy. Results show that relative changes in ctDNA positively correlate with tumor volume changes, with lower baseline ctDNA associated with better clinical responses. In some complete responders, ctDNA became undetectable before radiological confirmation, demonstrating both tumor reduction and ctDNA clearance. However, all patients eventually progressed. Parameter analysis revealed that treatment efficacy significantly impacts ctDNA shedding patterns, often causing characteristic peaks in ctDNA levels. These dynamics could serve as an early warning system for progression, potentially enabling more timely intervention. The model effectively characterizes patient-specific tumor and ctDNA dynamics. Results suggest alternative strategies, including chemotherapy, could optimize dosing regimens based on ctDNA patterns to improve responses and extend time to progression. This work establishes a foundation for integrating ctDNA surveillance into treatment monitoring for ASCC patients.
Biology-based modeling and topological data analysis (TDA) are powerful techniques for characterizing properties of tumors and vascular networks, but there has been limited effort to integrate these approaches to characterize in vivo tumor growth. Persistent homology offers a systematic approach to identify features such as connected components, loops, and voids across different scales in high-dimensional datasets. Likewise, biology-based models can be calibrated to longitudinal data to yield tumor-specific parameters describing tumor and vascular growth. In this study, we applied TDA and biology-based modeling to longitudinal MRI collected in nine animals with C6 glioma tumors. Animals were imaged up to seven times over a two week period to measure the cell density and blood volume fraction. We computed persistent homology of cubical complexes filtered by the ratio of blood volume fraction to normalized tumor cell density to characterise connected components, loops, and voids in the 3D data. We summarised the output in 15 topological features which quantify known biological properties of the data. For the biology-based modeling approach, a two-species reaction-diffusion model describing tumor growth and angiogenesis was calibrated to longitudinal data to estimate parameters describing tumor growth, invasion, angiogenesis, and vessel death. We then performed k-means clustering on a combined set of topological and biology-based modeling features yielding three clusters. Clusters 1 and 2 consisted of tumors that exhibited voids (necrosis), while Cluster 3 consisted of tumors without well-defined voids. Notably animals in Clusters 1 and 2 had a lower ratio of vascular proliferation to tumor proliferation than Cluster 3. This preliminary study indicates that there may be relationships between topological features and biology-based parameters. Further development of these methods could yield a framework to assign improved model parameters of tumor growth and response.
We begin with a system of ordinary differential equations for an immunological treatment of a primary tumor by neoantigen peptide vaccines. This system is coupled with a partial differential equation of metastasis that tracks the number of metastases per time and size. Vaccine dose is taken as a control in the primary tumor ordinary differential equation to slow tumor growth and the spread of metastatic tumors. An optimal control problem is formulated to design vaccine treatment.
Resistance to treatment, which comes from the heterogeneity of cell types within tumors, is a leading cause of poor treatment outcomes in cancer patients. Previous mathematical work modeling cancer over time has neither emphasized the relationship between cell heterogeneity and treatment resistance nor depicted heterogeneity with sufficient nuance. To respond to the need to depict a wide range of resistance levels, we develop a random differential equation model of tumor growth. Random differential equations are differential equations in which the parameters are random variables. In the inverse problem, we aim to recover the sensitivity to treatment as a probability mass function. This allows us to observe what proportions of cells exist at different sensitivity levels. After validating the method with synthetic data, we apply it to monoclonal and mixture cell population data of isogenic Ba/F3 murine cell lines to uncover each tumor's levels of sensitivity to treatment as a probability mass function. We emphasize the applications of this project by fitting the model to patient prostate cancer data to recover changes in treatment sensitivity over multiple treatment cycles.
Immunotherapies that target the host immune system to mount effective responses to cancer hold great promise. Overcoming patient- and organ-specific heterogeneity remain a significant challenge. In order to quantify individual patient responses to treatment, we fit a tumor-immune mathematical model to patient and site-specific dynamics of response to combination treatment (nivolumab + ipilimumab + entinostat) using tumor data (RECIST) coupled with immune markers measured by imaging mass cytometry. Bayesian parameter inference of the site-specific patient responses reveals that only immunosuppression parameters can explain response; parameters controlling cytotoxicity are not predictive. Via the fits of many tumors, we quantify the variability in tumor-immune dynamics across patients and tissues and reveal controllable parameter regimes. We go on to show that through posterior parameter sampling and simulation, we are able to use our model to extrapolate beyond the data and predict the probability of response in virtual metastatic tumors in patients for which we have no data at that site, thus overcoming the limitations of a small clinical trial to enable the analysis of a large virtual patient and virtual tumor cohort undergoing combination treatment.
CRISPR-based genome editing technologies have enabled massively-parallel genomic screens, such as DepMap – a Broad Consortium effort to catalog gene knockouts in cancer cell lines. These projects find that the growth effects of a mutation depend heavily on the background genotype of a cell. Evolutionary theory has studied the effects of background genotype on mutations for generations and has uncovered general patterns across the tree of life These patterns found in evolving populations have culminated in a ‘Geometric Model’ of adaptation that has successfully predicted the effects of novel combinations of mutations in yeast and E. coli. This model could in principle be applied to DepMap and other massively-parallel genomic screens to learn genotype to phenotype to fitness mappings and potentially predict the evolution of a population. Fitting this model to large-scale real data, however, is challenging because the model infers a latent (hidden) space of phenotypes with mathematical symmetries which confuse regression methods. Here, we present a methodology for fitting a Geometric Model of adaptation to large-scale genomic screens that eliminates rotational, translational, and permutation symmetries in the inferred phenotype space and successfully reconstructs genotype to phenotype to fitness mappings of Liver cancer cell line knockout data. Thus, making comprehensive quantitative models of genotype to phenotype to fitness mappings possible in a multitude of diseases, which in turn will allow us to infer phenotypic complexity and predict treatment response.
Mathematical oncology is an interdisciplinary research field in which mathematical modelling, analysis, and simulation are used to study cancer. In this work, we perform a bibliometric analysis to describe how mathematical oncology has changed over time. We quantitatively interrogate temporal trends in the field by analysing article metadata such as authors, publication dates, titles, article keywords, and abstracts. We specifically investigate if and how these trends have been shaped by paradigm-shifting research advances and world events. The data are collected from bibliographic databases such as Web of Science and Scopus, as well as the world's most prominent mathematical biology journals including: the Bulletin of Mathematical Biology, the Journal of Mathematical Biology, the Journal of Theoretical Biology, and Mathematical Biosciences. We show that, since the 1960's, mathematical oncology has become increasingly data-driven, international, and interdisciplinary.
Solid tumors such as pancreatic cancer are major causes of cancer-related deaths worldwide. Despite the availability of multiple treatment options such as radiotherapy or chemotherapy, long-term survival rates of patients with solid tumors remain low due to the development of treatment resistance and tumor recurrence. It has been experimentally observed that irradiation induces shifts in tumor growth kinetics, highlighting the need to unravel both short- and long-term cellular responses to irradiation. Computational models have been used to complement experimental studies by quantifying complex interactions between radiation, tumor biology, and treatment variables. While the common approach of employing the linear-quadratic model and its derivatives by computing the survival fraction is successful in describing short-term effects of radiation on a tumor, it is not suitable for capturing dynamic, persistent, long-term treatment effects. In this study, we developed a phenomenological differential equation-based model that integrates both immediate and delayed radiotherapy effects. A key feature of our model is the inclusion of probabilistic proliferation dynamics. We incorporate cancer cell proliferation rates as the determinant of radiosensitivity, aligning with the well-established hypothesis that highly proliferative cells are more radiosensitive than slower proliferating cells. By using these proliferation rates to determine the rate of cell death after irradiation, the model predicts a heterogeneous cancer cell killing rate, resulting in a variable fraction of surviving cells and a subsequent shift in the composition of the tumor. Thus, the model provides mechanistic insights into relapse dynamics and heterogeneous treatment responses. In the future, we want to extend our model by including immune cell dynamics to investigate the impact of radiation on the tumor microenvironment and the reciprocal interactions between cancer cells and the immune system.
Synthetic receptors are engineered proteins with potential applications to studying cell-cell interactions and cancer cell therapy. One promising research direction is engineering the Notch receptor, a transmembrane protein that can detect extracellular signals such as antigens or other ligands, and convert them to intracellular signals to activate expression of certain genes. Both the intracellular and extracellular domains can be engineered and replaced with alternative domains, creating the family of modified Notch receptors known as synthetic Notch (synNotch). SynNotch has the potential to improve chimeric antigen receptor (CAR) T-cell therapy by tuning binding affinity to a specific cancer antigen and minimizing off-target effects. We propose an ordinary differential equation model of synNotch receptor activity that has predictive value of how custom cell response behaviors can be programmed. The mathematical model couples activation dynamics on a fast timescale, characteristic of receptor-ligand interactions, and of slower downstream gene expression dynamics. Local and global sensitivity analyses indicate model parameters that yield the greatest variability in downstream results, indicating their potential to be engineered for different functions. Specifically, we find that ligand association and ligand-dependent activation have the greatest potential for modulating transcription factor release.
Radiotherapy planning typically assumes homogeneous efficacy across the tumor, which can lead to an overestimation of tumor cell death and control. We seek to quantify these errors by identifying the radiation dose boost required by non-homogeneous treatment efficacy models to yield the tumor response predicted by homogeneous response. We collected data from ten glioblastoma patients who received a total of 60 Gray (Gy) of radiation delivered in 30 fractions, concurrent chemotherapy, and were imaged prior to, and then at one-, three-, and five-months post-therapy. These data were used to inform a patient-specific, mechanically coupled reaction-diffusion model describing the spatiotemporal progression of tumor growth and response to therapy. The radiotherapy effect was modeled as an instantaneous decrease in the tumor volume fraction (ϕ) after each dose, with the surviving fraction (SF) defined as the ratio between the post- and pre-treatment ϕs. For each patient, the proliferation rates, diffusion coefficients, and SFs were calibrated to data up to five months post-therapy. We then simulated tumor growth using the calibrated model with (a) homogeneous SF, or heterogeneous SF, based on (b) vascularity or (c) cell density. We determined the patient-specific global dose boost required for models (b) and (c) to match the response predicted by model (a) one-month post-radiotherapy for SFs of 0.2-1.0 in increments of 0.05. For 0.55-0.95 SFs, model (b)’s predicted response required an additional 0.67 +/- 0.30 (mean +/- standard deviation) Gy per day, while model (c) only needed an additional 0.27 +/- 0.18 Gy per day across all patients. There was a statistically significant (P < 0.05, Wilcoxon rank sum test) difference between dose predictions for models (b) and (c) across all SFs. We developed an approach to calculate the radiation dose increases needed by non-homogeneous treatment efficacy models to match the tumor response predicted by a homogeneous model.
The analysis of prostate-specific antigen (PSA) dynamics after external beam radiotherapy is crucial for detecting prostate cancer recurrence. A significant increase in PSA post-radiotherapy often indicates biochemical relapse, although this evolution can be gradual and may take many years to manifest. Current clinical criteria for defining biochemical relapse rely on observation of population-based markers, using fixed thresholds to assess patient progression after a minimum value of PSA is reached. However, this approach does not account for individual tumor dynamics, which may delay recurrence detection and subsequent treatment. To overcome this limitation, we propose anticipating PSA increases using patient-specific forecasts obtained with a mechanistic model that describes post-radiotherapy tumor dynamics. This model utilizes longitudinal PSA measurements, which are routinely collected as part of standard-of-care management for prostate cancer before and after radiotherapy. By applying Bayesian calibration to the model using these data series, we can thus predict patient-specific PSA dynamics, accounting for the uncertainties in the model and data. Additionally, we can obtain the probabilistic distribution of key model-based biomarkers of biochemical relapse (e.g., surviving tumor cell proliferation rate, PSA nadir, and time to biochemical relapse), allowing for early identification of biochemically-relapsing patients. By leveraging our probabilistic formulation, we also introduce risk measures based on the distributions of these biomarkers to allow for a more accurate assessment of an individual’s risk of relapse (e.g., superquantiles). Finally, although validation in larger, more diverse cohorts is needed and extensions of the model could be implemented, this approach has the potential to improve clinical decision-making by personalizing the monitoring of prostate cancer patients after radiotherapy and anticipating disease progression to advanced stages.
Non-muscle-invasive bladder cancer (NMIBC) is one of the most prevalent oncological diseases worldwide, originating in the bladder epithelium, known as the urothelium. Mitomycin C (MMC) chemotherapy is a widely used treatment that reduces recurrence rates and prolongs progression-free survival. However, its full mechanism of action in BC and its immune-related effects, which are crucial for the formulation of an ideal regimen of MMC, remains to be elucidated. This work integrates systems immunology principles with temporal ordinary differential equations (ODEs) to provide a test bed for the theoretical investigation of immune system dynamics during disease progression and chemotherapy administration. We first identify distinct tumor and immune populations and formulate their specific interactions based on biological research. After, we simulate hypothetical BC cases to illustrate the complex dynamics of specialized cell types that bridge the innate and adaptive immune responses. [1] Yosef, M., & Bunimovich-Mendrazitsky, S. (2024). Mathematical model of MMC chemotherapy for non-invasive bladder cancer treatment. Frontiers in Oncology, 14. https://doi.org/10.3389/fonc.2024.1352065 [2] Bunimovich-Mendrazitsky, S., Pisarev, V., & Kashdan, E. (2015). Modeling and simulation of a low-grade urinary bladder carcinoma. Computers in Biology and Medicine, 58, 118–129. https://doi.org/10.1016/j.compbiomed.2014.12.022
This study uses a vertex-based computational model based on CHASTE to investigate interclonal cooperation within cell populations. We analyze interactions between mutated cell types with distinct hyperproliferative and invasive properties. Modeling the cooperative dynamics between these clones shows how cellular heterogeneity affects tissue growth and invasion, enhancing the understanding of tumor progression and therapeutic target predictions.
Background: Patient-reported outcomes (PROs) are defined by the FDA as “any report of the status of a patient's health condition that comes directly from the patient, without interpretation of the patient's response by a clinician or anyone else.” PROs are promising investigative biomarkers for cancer response and progression due to their being non-invasive and low-cost as well as their ease and frequency to be administered to patients. In this study, we seek to elucidate the relationship between tumor size (TS) and PRO dynamics and leverage this information to infer patient-specific tumor growth inhibition (TGI) parameters. Methods: We developed a joint model describing longitudinal PROs using a Markov chain model with TS as a time-dependent covariate in the transition rate matrix. We then modeled TS using the Claret TGI model. We then simulated 10 in silico patients. To test how the informativeness of the simulated PROs, we performed Bayesian inference on the probability distribution of TGI parameters given: 1) no data as input); 2) only simulated TS data; 3) only simulated PRO data; and 4) both TS and PRO data. Results: Considering PRO data with or without TS data increased precision (1 / standard deviation) of patient-specific TGI parameters by a factor of 1.9 (1.4–2.6) and 1.1 (0.7–1.8), respectively. By contrast, considering PRO data with or without TS data only marginally improved accuracy (1 / root mean squared error) of patient-specific TGI parameters by a factor of 1.4 (0.4–4.5) and 1.3 (0.5–3.5). Conclusion: These results suggest that PROs have the promise to be leveraged as a minimally invasive and inexpensive biomarker to inform patient-specific TGI parameters. Future research directions include exploring the effect of sampling frequency, number of patients, number of PROs, and model misspecification, as well as application on a clinical dataset of 63 NSCLC patients treated with immune checkpoint inhibitors.
In this work, we address the treatment of glioblastoma (GBM), a difficult-to-treat brain cancer. Upon discovery of GBM, patients are frequently treated with both surgery and chemoradiotherapy but still suffer from eventual recurrence due to unresectable microscopic disease that evades adjuvant therapy. The disease is characterized by an immunosuppressive tumor microenvironment. Immunotherapies, including immune checkpoint inhibitors, have been trialed in GBM but have so far failed to improve otherwise bleak outcomes for GBM patients. One possible way GBM tumors sustain this immune suppression, even in the face of ICI, is through recruitment and sustaining a population of myeloid-derived suppressor cells (MDSCs). These potently immunosuppressive cells decrease T-cell activity through a range of mechanisms including direct cell-cell interactions, secretion of T-cell inhibitors, and alteration of the metabolic environment. Furthermore, radiotherapy (RT), which can paradoxically lead to both immunosuppression and immune stimulation, represents an underexplored option to possibly tip the tumor in favor of immune stimulation. To understand the interplay between MDSCs, effector cells, and RT, we developed a dynamic, computational model of the tumor immune microenvironment of glioblastoma.
We present a model for T-cell re-activation under the action of therapeutic compounds targeting the PD-1/PD-L1 interaction. The effector T cells in the assay express Luciferase upon TCR-mediated activation, which is diminished by the presence of PD-1/PD-L1 interaction provided by co-cultured artificial antigen-presenting cells. Upon PD-L1 blockade with tested compounds added for 6 hours at different dose concentrations, the activation of effector T cells is restored, reflected by increased luminescence signals. The resulting dose response curve is non-monotone due to the existence of two separate and antagonizing effects - specific activation of T cells and unspecific toxicity, observed separately, but also overlapping at a certain range of the compound concentration. A mathematical model is used to estimate the concentration for maximum level of activation and the EC50 concentration. The model presents a mechanism for the temporal change in the activity of the Firefly Luciferase over the course of the experiment as measured by the strength of the luminescent signal based on a system of ODEs. Parameters of the model are estimated from experimental data and are used to estimate the drug concentration corresponding to the maximum level of T-cell activation and the EC50 concentration. This is joint work with Lukasz Skalniak (Department of Organic Chemistry, Faculty of Chemistry, Jagiellonian University, Krakow, Poland). This work was supported by grant DP-05-KOST-13 from FNI, Bulgarian Ministry of Education and Science and by grant 2021/42/E/NZ7/00422 from the National Science Centre, Poland.
Circulating tumor DNA (ctDNA) has emerged as a promising biomarker for monitoring cancer progression and treatment response in real time. In anal squamous cell carcinoma (ASCC), where 80-90% of cases are linked to human papillomavirus (HPV), ctDNA demonstrates high sensitivity in tracking disease dynamics, often detecting progression earlier than imaging while enabling frequent assessment and correlating with tumor burden. Our study examined how patient-specific modeling of ctDNA dynamics can predict time to progression in HPV-associated ASCC. We analyzed longitudinal data from 32 ASCC patients receiving immunotherapy every 3 weeks for up to 2 years, exploring correlations between tumor volume and ctDNA levels. We developed a mathematical model calibrated to patient-specific tumor volume and ctDNA dynamics during immunotherapy. Results show that relative changes in ctDNA positively correlate with tumor volume changes, with lower baseline ctDNA associated with better clinical responses. In some complete responders, ctDNA became undetectable before radiological confirmation, demonstrating both tumor reduction and ctDNA clearance. However, all patients eventually progressed. Parameter analysis revealed that treatment efficacy significantly impacts ctDNA shedding patterns, often causing characteristic peaks in ctDNA levels. These dynamics could serve as an early warning system for progression, potentially enabling more timely intervention. The model effectively characterizes patient-specific tumor and ctDNA dynamics. Results suggest alternative strategies, including chemotherapy, could optimize dosing regimens based on ctDNA patterns to improve responses and extend time to progression. This work establishes a foundation for integrating ctDNA surveillance into treatment monitoring for ASCC patients.
Biology-based modeling and topological data analysis (TDA) are powerful techniques for characterizing properties of tumors and vascular networks, but there has been limited effort to integrate these approaches to characterize in vivo tumor growth. Persistent homology offers a systematic approach to identify features such as connected components, loops, and voids across different scales in high-dimensional datasets. Likewise, biology-based models can be calibrated to longitudinal data to yield tumor-specific parameters describing tumor and vascular growth. In this study, we applied TDA and biology-based modeling to longitudinal MRI collected in nine animals with C6 glioma tumors. Animals were imaged up to seven times over a two week period to measure the cell density and blood volume fraction. We computed persistent homology of cubical complexes filtered by the ratio of blood volume fraction to normalized tumor cell density to characterise connected components, loops, and voids in the 3D data. We summarised the output in 15 topological features which quantify known biological properties of the data. For the biology-based modeling approach, a two-species reaction-diffusion model describing tumor growth and angiogenesis was calibrated to longitudinal data to estimate parameters describing tumor growth, invasion, angiogenesis, and vessel death. We then performed k-means clustering on a combined set of topological and biology-based modeling features yielding three clusters. Clusters 1 and 2 consisted of tumors that exhibited voids (necrosis), while Cluster 3 consisted of tumors without well-defined voids. Notably animals in Clusters 1 and 2 had a lower ratio of vascular proliferation to tumor proliferation than Cluster 3. This preliminary study indicates that there may be relationships between topological features and biology-based parameters. Further development of these methods could yield a framework to assign improved model parameters of tumor growth and response.
We begin with a system of ordinary differential equations for an immunological treatment of a primary tumor by neoantigen peptide vaccines. This system is coupled with a partial differential equation of metastasis that tracks the number of metastases per time and size. Vaccine dose is taken as a control in the primary tumor ordinary differential equation to slow tumor growth and the spread of metastatic tumors. An optimal control problem is formulated to design vaccine treatment.
Resistance to treatment, which comes from the heterogeneity of cell types within tumors, is a leading cause of poor treatment outcomes in cancer patients. Previous mathematical work modeling cancer over time has neither emphasized the relationship between cell heterogeneity and treatment resistance nor depicted heterogeneity with sufficient nuance. To respond to the need to depict a wide range of resistance levels, we develop a random differential equation model of tumor growth. Random differential equations are differential equations in which the parameters are random variables. In the inverse problem, we aim to recover the sensitivity to treatment as a probability mass function. This allows us to observe what proportions of cells exist at different sensitivity levels. After validating the method with synthetic data, we apply it to monoclonal and mixture cell population data of isogenic Ba/F3 murine cell lines to uncover each tumor's levels of sensitivity to treatment as a probability mass function. We emphasize the applications of this project by fitting the model to patient prostate cancer data to recover changes in treatment sensitivity over multiple treatment cycles.
CRISPR-based genome editing technologies have enabled massively-parallel genomic screens, such as DepMap – a Broad Consortium effort to catalog gene knockouts in cancer cell lines. These projects find that the growth effects of a mutation depend heavily on the background genotype of a cell. Evolutionary theory has studied the effects of background genotype on mutations for generations and has uncovered general patterns across the tree of life These patterns found in evolving populations have culminated in a ‘Geometric Model’ of adaptation that has successfully predicted the effects of novel combinations of mutations in yeast and E. coli. This model could in principle be applied to DepMap and other massively-parallel genomic screens to learn genotype to phenotype to fitness mappings and potentially predict the evolution of a population. Fitting this model to large-scale real data, however, is challenging because the model infers a latent (hidden) space of phenotypes with mathematical symmetries which confuse regression methods. Here, we present a methodology for fitting a Geometric Model of adaptation to large-scale genomic screens that eliminates rotational, translational, and permutation symmetries in the inferred phenotype space and successfully reconstructs genotype to phenotype to fitness mappings of Liver cancer cell line knockout data. Thus, making comprehensive quantitative models of genotype to phenotype to fitness mappings possible in a multitude of diseases, which in turn will allow us to infer phenotypic complexity and predict treatment response.
Mathematical oncology is an interdisciplinary research field in which mathematical modelling, analysis, and simulation are used to study cancer. In this work, we perform a bibliometric analysis to describe how mathematical oncology has changed over time. We quantitatively interrogate temporal trends in the field by analysing article metadata such as authors, publication dates, titles, article keywords, and abstracts. We specifically investigate if and how these trends have been shaped by paradigm-shifting research advances and world events. The data are collected from bibliographic databases such as Web of Science and Scopus, as well as the world's most prominent mathematical biology journals including: the Bulletin of Mathematical Biology, the Journal of Mathematical Biology, the Journal of Theoretical Biology, and Mathematical Biosciences. We show that, since the 1960's, mathematical oncology has become increasingly data-driven, international, and interdisciplinary.
Solid tumors such as pancreatic cancer are major causes of cancer-related deaths worldwide. Despite the availability of multiple treatment options such as radiotherapy or chemotherapy, long-term survival rates of patients with solid tumors remain low due to the development of treatment resistance and tumor recurrence. It has been experimentally observed that irradiation induces shifts in tumor growth kinetics, highlighting the need to unravel both short- and long-term cellular responses to irradiation. Computational models have been used to complement experimental studies by quantifying complex interactions between radiation, tumor biology, and treatment variables. While the common approach of employing the linear-quadratic model and its derivatives by computing the survival fraction is successful in describing short-term effects of radiation on a tumor, it is not suitable for capturing dynamic, persistent, long-term treatment effects. In this study, we developed a phenomenological differential equation-based model that integrates both immediate and delayed radiotherapy effects. A key feature of our model is the inclusion of probabilistic proliferation dynamics. We incorporate cancer cell proliferation rates as the determinant of radiosensitivity, aligning with the well-established hypothesis that highly proliferative cells are more radiosensitive than slower proliferating cells. By using these proliferation rates to determine the rate of cell death after irradiation, the model predicts a heterogeneous cancer cell killing rate, resulting in a variable fraction of surviving cells and a subsequent shift in the composition of the tumor. Thus, the model provides mechanistic insights into relapse dynamics and heterogeneous treatment responses. In the future, we want to extend our model by including immune cell dynamics to investigate the impact of radiation on the tumor microenvironment and the reciprocal interactions between cancer cells and the immune system.
Many foraging animals, including bees, develop near-optimal movement patterns based on memory. While models have simulated how bees establish deterministic traplines, formal mathematical proofs of their behavior remain scarce. We address this gap by adapting and simplifying the Dubois et al. (2024) model to enable mathematical analysis. We prove that simulated bees will always eventually converge to a single deterministic route. Additionally, we propose conjectures about the distribution of routes to which simulated bees may converge. Future work could explore inference methods for learned behavior based on this model. These findings have implications beyond biology, providing insights into reinforced random walks and reinforcement learning.
In this work, we study the applicability of the Complete Electrode Model (CEM) to Electroencephalography (EEG). EEG is a non-invasive imaging technique where it aims to localize cerebral sources generating the measured EEG signals. An existence and uniqueness result of this model is proved. Preliminary numerical implementation are also done.
Chemical Reaction Networks provide a fundamental framework for modeling the stochastic dynamics of biochemical systems, where molecular species evolve through discrete and random noise reaction events. Parameter inference in Chemical Reaction Networks is a central prob- lem in systems biology, but traditional methods such as maximum likelihood estimation are often intractable due to computational complexity and the lack of continuous-time data. In this study, we introduce a statistically grounded and computationally efficient estimator for reac- tion rate parameters using high-frequency discrete-time observations. Modeling the system as a Continuous-Time Markov Chain, our method handles general kinetics, including non-mass- action and higher-order reactions. Validation on synthetic and experimental datasets demon- strates its accuracy and robustness. This approach offers a simple and reliable framework for parameter inference in complex stochastic systems.
Many foraging animals, including bees, develop near-optimal movement patterns based on memory. While models have simulated how bees establish deterministic traplines, formal mathematical proofs of their behavior remain scarce. We address this gap by adapting and simplifying the Dubois et al. (2024) model to enable mathematical analysis. We prove that simulated bees will always eventually converge to a single deterministic route. Additionally, we propose conjectures about the distribution of routes to which simulated bees may converge. Future work could explore inference methods for learned behavior based on this model. These findings have implications beyond biology, providing insights into reinforced random walks and reinforcement learning.
In this work, we study the applicability of the Complete Electrode Model (CEM) to Electroencephalography (EEG). EEG is a non-invasive imaging technique where it aims to localize cerebral sources generating the measured EEG signals. An existence and uniqueness result of this model is proved. Preliminary numerical implementation are also done.
SMB2025