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.
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.
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.
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).
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.
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.
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.
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.
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.
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.
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.
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
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.
SMB2025