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

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

The retinal microcirculation supplies blood and oxygen to the cells responsible for vision, and vascular impairments – including compromised flow regulation, reduced capillary density, and elevated intraocular pressure – are involved in the progression of eye diseases such as glaucoma. Here, an established theoretical model of a retinal microvascular network will be presented and extended to investigate the effects of these impairments on retinal blood flow and oxygenation as intraluminal pressure is varied. A heterogeneous description of the arterioles based on confocal microscopy images is combined with a compartmental representation of the downstream capillaries and venules. A Green’s function method is used to simulate oxygen transport in the arterioles, and a Krogh cylinder model is used in the capillary and venular compartments. Acute blood flow autoregulation is simulated in response to changes in pressure, shear stress, and metabolism. The model predicts that impaired flow regulation mechanisms, decreased capillary density, and increased intraocular pressure all cause a loss in the autoregulation plateau over the baseline range of intraluminal pressures (meaning that blood flow is not maintained constant over those pressures), leading to a corresponding decrease in oxygenation in that range. Small impairments in capillary density or intraocular pressure are predicted to mostly be offset by functional flow regulation; however, larger changes and/or combinations of vascular impairments lead to a significant decrease in oxygenation. Clinically, since poor retinal tissue oxygenation could lead to vision loss in advanced glaucoma, model results suggest early identification of vascular changes to prevent these impairments from progressing.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Cells migrate in response to gradients in extra-cellular chemical signals in a process known as chemotaxis. Recent experiments on the model microorganism Dictyostelium discoideum have shown that dense aggregates of cells collectively undergoing chemotaxis exhibit emergent fluid-like properties such as viscosity and surface tension. In this work, we use simulations to explain how active interactions between cells give rise to these emergent phenomena. We propose an agent-based model for intermittent cell-cell attachments and show that it gives rise to emergent fluid-like behavior for an aggregate of cells. We generalize this model to include cell-surface attachments, and show that surface-associated aggregates display properties similar to a liquid droplet resting on a surface. Furthermore, we study the situation where cells self-generate and respond to a chemical gradient by consuming an externally supplied chemoattractant. Our simulations reveal how individual cells move inside the swarm as the cells move as a collective. Finally, we predict some of the key cellular processes that are responsible for this collective behavior, and provide hypotheses to be tested in future experimental studies.

Chemical signaling, or quorum sensing (QS), promotes collective behaviour in bacteria, from biofilm formation to swarming. By coupling QS systems with genes that control motility, bacteria can be engineered to generate tunable spatio-temporal patterns in vitro. However, it is not well-understood in general how the type of gene-regulatory network affects emergent population-level patterning. In this talk, we investigate the effect of the gene-regulatory network on emergent patterning in a population of motile bacteria that interact via QS. By formally upscaling a cell-level model in a biologically relevant scaling regime, we derive a continuum model that explicitly accounts for genetic regulation of motility and signal production through chemical structuring. Using a WKBJ-like framework, we derive criteria for the onset of two types of emergent patterning for a canonical QS circuit. We also uncover a new route to the well-known phenomenon of motility-induced phase separation (MIPS) through genetic regulation of tumbling frequency. Lastly, we discuss generalisations of our WKBJ-like analysis to more complex gene-regulatory networks that exhibit bistability.

Stem cells maintain and repair our tissues, but not all stem cells are identical. As organisms age, distinct stem cell 'clones' can begin to dominate the cell population. While this behavior has been observed across multiple species and organs, the mechanisms and consequences of stem cell clonality are still poorly understood. We have developed a novel experimental approach using a CRISPR-Cas9 system to uniquely “barcode” spermatogonial stem cell clones in the testes of male zebrafish. Once these fish reach sexual maturity, we sample sperm each month to determine the contribution of each stem cell clone to the sperm pool over time. The observed clonal dynamics may be driven by factors such as genetic drift, selection, or sampling error. We hypothesize that a small number of clones are under positive selection, resulting in their eventual dominance in the sperm pool. To bridge the gap between theory and data, we have developed stochastic models of stem cell dynamics in the testis. These models are formulated as hidden Markov models that describe rules for the division and differentiation of stem cells within the testis. We first evaluate our ability to estimate model parameters on simulated data. Then, we apply our model to the experimental data to quantify evidence for genetic drift and selection. Our models provide insight into how individual stem cell behavior can lead to population-level clonality.

Osteoporosis is a disease that affects more than 200 million people all over the world. Although its underlying mechanisms are gradually being revealed, effective treatments or preventive measures have not been established yet. This study focused on age-related osteoporosis by measuring bone mass and bone metabolism markers in mice from 4 to 52 weeks of age. We developed mathematical models describing bone metabolism and analyzed experimental data. From the result of data analysis, we quantitatively elucidated the mechanisms of bone loss. Furthermore, we conducted treatment intervention simulations by changing parameter values in mathematical models to identify effective bone metabolism pathways for increasing bone mass and new potential therapeutic strategies.

Retinal detachment occurs when the neurosensory retina separates from the retinal pigment epithelium (RPE), disrupting the nutrient supply to photoreceptor cells. There are three types of retinal detachment: exudative (ERD), tractional (TRD), and rhegmatogenous (RRD), with RRD being the most common. RRD develops when a retinal tear or hole allows vitreous humor to enter the subretinal space, causing the neurosensory retina to detach from the RPE. If left untreated, this condition can lead to irreversible vision loss. Although ophthalmological tools can detect RRD, its rate of progression—particularly due to continuous eye movement—remains poorly understood. This study develops a fluid-structure interaction model to examine how various factors, including retinal thickness, elasticity, adhesion strength between the retina and RPE, vitreous humor density and viscosity, and eye rotation speed, influence detachment progression. By quantifying detachment rates under different conditions, this research aims to enhance our understanding of RRD dynamics and refine estimates of effective treatment timelines to prevent permanent visual impairment. Student: William Ebo Annan Advisors: Prof. Diana White & Prof. Emmanuel O.A. Asamani

Genetically identical cells can adopt distinct, stable states, playing a crucial role in development and tissue organisation. This talk explores the mechanisms driving cell fate decisions, focusing on the interplay between gene regulatory networks and cell-to-cell communication. Using spatial stochastic models that capture fine-scale regulatory dynamics, we demonstrate how feedback loops and paracrine signalling function as switch-like controllers of cell fate, enabling transitions from transient to stable states. We derive mathematical expressions predicting the threshold signalling strength required to trigger phase transitions and establish a fundamental limit on the spatial spread of phenotypic regions. Specifically, we show that the mean region size scales proportionally to the cubic root of the signalling strength, implying that large, stable domains are prohibitively costly to maintain. This trade-off between robustness and signalling precision highlights the constraints organisms must navigate during development to maintain spatial organisation. Our findings provide key insights into the principles governing multicellular patterning and the regulation of tissue structure.

Multiscale problems are prevalent in many real world scenarios, especially in biology, where the behaviour of a single microorganism can have considerable impact over lengthscales much larger than its own. In this work we consider the effect of the membrane microstructure of a bacterial cell on the behaviour of concentration profiles of relevant molecules on bacterial and bacterial colony lengthscales. Transport through the outer membrane of gram-negative bacteria is restricted to specific channels and non-specific porins. These provide a size-restricted passageway for small molecules through an otherwise impermeable membrane. The effects of these channels are important, for example, antibiotics must cross the outer membrane in order to effectively target gram-negative bacteria, and quorum sensing molecules must cross the membrane to allow bacterial colonies to coordinate mass phenotypic changes such as the production of virulence factors. In mathematical models this limiting transport mechanism across the membrane is often represented via phenomenological constitutive boundary conditions. In this work, we systematically derive the correct effective boundary conditions to impose across a bacterial membrane in terms of physical channel and porin properties. We use a hybrid mathematical approach, combining multiscale methodology such as asymptotic homogenisation and boundary layer theory with numerical simulations. More broadly, because we consider a generic membrane geometry and do not impose a specific outer problem, the results that we derive have a wide scope of potential applications beyond bacterial membranes, for example, to model water vapour or heat loss through fabrics, or mass transfer through surface coatings in chemical engineering.

Many genes are expressed in bursts of transcription, associated with alternating active and inactive promoter states. Such transcriptional bursting is characterised by the burst frequency and burst size, which describe how often a burst occurs and how many transcripts are produced per burst. These two burst parameters offer a simple, intuitive and practical quantitative description of bursty gene expression dynamics. However, a transcriptome-wide picture of how the burst frequency and size are modulated due to gene replication and other cell-cycle dependent factors remains missing. To address this, we fit mechanistic models of gene expression to mRNA count data for thousands of mouse genes, obtained by sequencing of single cells whose cell-cycle position has been inferred previously. Although we observe substantial heterogeneity in transcriptional regulation, we find that upon DNA replication, the genome-wide median burst frequency approximately halves, while the median burst size remains mostly unchanged, thus shedding light on the effect of gene dosage compensation. We show that to accurately estimate the bursting kinetics from sequencing data, mechanistic models must explicitly account for gene copy number variation and extrinsic noise due to factors varying across the cell cycle, whereas correcting for technical noise due to imperfect mRNA capture is less critical.

Pre-mRNA often contains introns, which are non-coding sequences that need to be cut out or spliced before translation occurs. The spliceosome, an essential catalyst composed of several proteins with specific sequence affinities, is required for this process. Very long introns must be removed in pieces, a process known as recursive splicing. The experimental literature on the time it takes for the splicing process to occur is inconsistent. Splicing was traditionally believed to be a slow process that could take anywhere from one to tens of minutes per splicing event. However, recent reports suggest that some splicing events occur within a few tens of seconds. We developed the chemical master equation corresponding to the biochemical mechanism of splicing, allowing us to derive the system’s probability distribution, and perform a stability analysis on two conditions based on an unknown association constant parameter associated with the binding step of the scaffolding complex. We also concluded that the distribution of splice times for a single event ranges from a few tens of seconds to a few tens of minutes. Through sensitivity analyses, we have found that the mean splicing time and distribution are almost entirely dependent on the rate at which the spliceosome is activated in the assembly process—i.e. when the U1 and U4 splicing factors dissociate—which confirms that this is the rate limiting step in the catalytic process. Finally, we have examined the distributions of recursive splicing up to six events, and derived analytic solutions for these recursive splicing events in the case where the scaffolding complex strongly binds to the pre-mRNA complex (the condition thought to favor recursive binding), thus providing a model that can be fit to experimental data to in order to evaluate the number of recursive splicing events occurring.

This work examines the interplay between age-structure and migration dynamics in collective cell behaviour. We focus on the integration of cell cycle dynamics with spatial migration, particularly examining the 'go-or-grow' hypothesis in the context of age-dependent processes. Our framework extends classical travelling wave theory to account for the age structure of cell populations, offering new insights into how cell cycle phases influence moving fronts and invasion dynamics. We analyse wave speed characteristics and front dynamics in age-structured systems, addressing a significant gap in current mathematical biology literature. The research provides a novel theoretical foundation for understanding how cell-cycle dependent proliferation and migration behaviours contribute to collective cell dynamics.

Cells respond to their local environment through mechanotransduction, converting mechanical signals into a biological response (e.g. cell growth, proliferation or differentiation). The cell cytoskeleton, particularly actomyosin stress fibres (SFs), and focal adhesions (FAs), which bind the cytoskeleton to the extra-cellular matrix (ECM), are central to this process, activating intracellular signalling cascades in response to deformation. We present a novel two-dimensional bio-chemo-mechanical model to describe the development of these structures, coupled through a positive feedback loop, and the associated cell deformation. Building on our previous one-dimensional approach, we similarly employ reaction-diffusion-advection equations to describe the evolution of key scaffolding and signalling proteins, and connect their concentrations to a viscoelastic description of the cell cytoplasm, ECM and adhesions. Further, we now incorporate other key mechanotransducing structures including the stiff cell nucleus, and plasma and cortical membranes. Working in an axisymmetric framework, we employ this model to explain how, dependent upon the mechanical properties of the surrounding ECM, non-uniform patterns of cell striation develop, leading to FA and SF localisation at the cell periphery. Moreover, a linear stability analysis reveals the stability of the axisymmetric configuration to various normal modes of deformation. By identifying non-axisymmetric modes with positive growth rates our model demonstrates a possible mechanism for self-driven surface patterning of cells in vitro.

In a number of bacteria, nitric oxide (NO) is converted to nitrate by an enzyme called Hmp. In emph{Streptomyces coelicolor}, synthesis of Hmp is in turn controlled by an iron-sulfur protein called NsrR. NsrR represses the transcription of two copies of the emph{hmp} gene in the emph{S. coelicolor} genome, but reaction of NsrR's iron-sulfur cluster with NO causes NsrR to dissociate from the emph{hmp} promoter, thus allowing Hmp to be expressed. While this is a straightforward control mechanism, NsrR is a dimer, and the iron-sulfur cluster in each monomer of NsrR can react with NO several times. Eventually, a repair system restores the NO-damaged iron-sulfur clusters of the dimers. But given that a single reaction with NO is sufficient to cause the NsrR dimer to dissociate from the emph{hmp} promoter, do we need to model the complex chemistry of the dimer, or is a highly simplified model that considers a single NsrR unit and its iron-sulfur cluster sufficient to capture the dynamics of this control system?

Across the spectrum of cell types, the concentration of calcium controls a wide array of cellular functions. These calcium signals, usually in the form of periodic oscillations, play a paramount role in correct cellular activity. T cells are fundamental to the correct behaviour of the immune system. These cells have recently been shown to exhibit two competing oscillatory mechanisms, depending on the influx of extracellular Ca2+. Ca2+ influx is controlled by two molecules, STIM1 and STIM2. When both STIMs are present, T cells showcase sinusoidal Ca2+ oscillations on a raised baseline, but when one of them is absent the nature of the oscillation changes to a mix of Ca2+ spikes and bursting periods. In this talk, we will present an ODE that attempts to explain how these two molecules control the nature of these oscillations in T cells

The prolactin receptor (PRLR) is a single-pass transmembrane receptor driving pancreatic β-cell proliferation via JAK/STAT signaling activation. This signal transduction pathway enables insulin-secreting β-cells to adapt to metabolic stress; however, the precise mechanisms underlying the pathway’s proliferative effect remain ill-defined. Here we implement a pipeline that uses live-cell fluorescence imaging, reconstitution approaches, and fluorescence correlation spectroscopy (FCS) to inform a mathematical model of PRLR signaling in β-cells and build a quantitative, mechanistic understanding of the signaling network. PRLR signaling is dynamic, involving changes in the spatial organization of signaling molecules. We have observed PRLR undergoing rapid internalization, a behavior that has been shown and modeled in other signaling pathways but has not been considered in a mathematical model of PRLR signaling. Such a model is useful for predicting strategies to modulate β-cell function. PRLR internalization is observed in both our minimal engineered PRLR expression system and in native pancreatic tissue, while FCS and chemigenetic labeling with SNAP-tag confirm the presence of a low concentration plasma membrane pool of PRLR. Our imaging data are used to integrate PRLR trafficking dynamics into an ordinary differential equation (ODE) model of PRLR signaling. We employ the ODE model to test hypotheses targeting how the spatial heterogeneity of PRLR signaling dynamics affects downstream signaling outcomes. Our data underscore the versatility of building a generalizable modeling-imaging framework to quantitatively understand signal transduction in and beyond β-cells.

Collective cell migration plays a crucial role in numerous biological processes, including cancer growth, wound healing, and the immune response. Often, the migrating population consists of cells with various different phenotypes. This study derives a general mathematical framework for modelling cell migration into the micro-environment, which is coarse-grained from an underlying individual-based model that captures some of the dynamics of cell migration that are influenced by the phenotype of the cell, such as: random movement, proliferation, phenotypic transitions, and interactions with the external environment. The resulting model provides a continuum, macroscopic description of cell invasion, which represents the phenotype of the cell as a continuous variable and is much more amenable to simulation and analysis than its individual-based counterpart when considering a large number of phenotypes. The results highlight how phenotypic structuring impacts the spatial and temporal dynamics of cell populations, demonstrating that different environmental pressures and phenotypic transition mechanisms significantly influence invasion patterns.

Cells respond to their local environment through mechanotransduction, converting mechanical signals into a biological response (e.g. cell growth, proliferation or differentiation). The cell cytoskeleton, particularly actomyosin stress fibres (SFs), and focal adhesions (FAs), which bind the cytoskeleton to the extra-cellular matrix (ECM), are central to this process, activating intracellular signalling cascades in response to deformation. We present a novel two-dimensional bio-chemo-mechanical model to describe the development of these structures, coupled through a positive feedback loop, and the associated cell deformation. Building on our previous one-dimensional approach, we similarly employ reaction-diffusion-advection equations to describe the evolution of key scaffolding and signalling proteins, and connect their concentrations to a viscoelastic description of the cell cytoplasm, ECM and adhesions. Further, we now incorporate other key mechanotransducing structures including the stiff cell nucleus, and plasma and cortical membranes. Working in an axisymmetric framework, we employ this model to explain how, dependent upon the mechanical properties of the surrounding ECM, non-uniform patterns of cell striation develop, leading to FA and SF localisation at the cell periphery. Moreover, a linear stability analysis reveals the stability of the axisymmetric configuration to various normal modes of deformation. By identifying non-axisymmetric modes with positive growth rates our model demonstrates a possible mechanism for self-driven surface patterning of cells in vitro.

In a number of bacteria, nitric oxide (NO) is converted to nitrate by an enzyme called Hmp. In emph{Streptomyces coelicolor}, synthesis of Hmp is in turn controlled by an iron-sulfur protein called NsrR. NsrR represses the transcription of two copies of the emph{hmp} gene in the emph{S. coelicolor} genome, but reaction of NsrR's iron-sulfur cluster with NO causes NsrR to dissociate from the emph{hmp} promoter, thus allowing Hmp to be expressed. While this is a straightforward control mechanism, NsrR is a dimer, and the iron-sulfur cluster in each monomer of NsrR can react with NO several times. Eventually, a repair system restores the NO-damaged iron-sulfur clusters of the dimers. But given that a single reaction with NO is sufficient to cause the NsrR dimer to dissociate from the emph{hmp} promoter, do we need to model the complex chemistry of the dimer, or is a highly simplified model that considers a single NsrR unit and its iron-sulfur cluster sufficient to capture the dynamics of this control system?

Across the spectrum of cell types, the concentration of calcium controls a wide array of cellular functions. These calcium signals, usually in the form of periodic oscillations, play a paramount role in correct cellular activity. T cells are fundamental to the correct behaviour of the immune system. These cells have recently been shown to exhibit two competing oscillatory mechanisms, depending on the influx of extracellular Ca2+. Ca2+ influx is controlled by two molecules, STIM1 and STIM2. When both STIMs are present, T cells showcase sinusoidal Ca2+ oscillations on a raised baseline, but when one of them is absent the nature of the oscillation changes to a mix of Ca2+ spikes and bursting periods. In this talk, we will present an ODE that attempts to explain how these two molecules control the nature of these oscillations in T cells

The prolactin receptor (PRLR) is a single-pass transmembrane receptor driving pancreatic β-cell proliferation via JAK/STAT signaling activation. This signal transduction pathway enables insulin-secreting β-cells to adapt to metabolic stress; however, the precise mechanisms underlying the pathway’s proliferative effect remain ill-defined. Here we implement a pipeline that uses live-cell fluorescence imaging, reconstitution approaches, and fluorescence correlation spectroscopy (FCS) to inform a mathematical model of PRLR signaling in β-cells and build a quantitative, mechanistic understanding of the signaling network. PRLR signaling is dynamic, involving changes in the spatial organization of signaling molecules. We have observed PRLR undergoing rapid internalization, a behavior that has been shown and modeled in other signaling pathways but has not been considered in a mathematical model of PRLR signaling. Such a model is useful for predicting strategies to modulate β-cell function. PRLR internalization is observed in both our minimal engineered PRLR expression system and in native pancreatic tissue, while FCS and chemigenetic labeling with SNAP-tag confirm the presence of a low concentration plasma membrane pool of PRLR. Our imaging data are used to integrate PRLR trafficking dynamics into an ordinary differential equation (ODE) model of PRLR signaling. We employ the ODE model to test hypotheses targeting how the spatial heterogeneity of PRLR signaling dynamics affects downstream signaling outcomes. Our data underscore the versatility of building a generalizable modeling-imaging framework to quantitatively understand signal transduction in and beyond β-cells.

Collective cell migration plays a crucial role in numerous biological processes, including cancer growth, wound healing, and the immune response. Often, the migrating population consists of cells with various different phenotypes. This study derives a general mathematical framework for modelling cell migration into the micro-environment, which is coarse-grained from an underlying individual-based model that captures some of the dynamics of cell migration that are influenced by the phenotype of the cell, such as: random movement, proliferation, phenotypic transitions, and interactions with the external environment. The resulting model provides a continuum, macroscopic description of cell invasion, which represents the phenotype of the cell as a continuous variable and is much more amenable to simulation and analysis than its individual-based counterpart when considering a large number of phenotypes. The results highlight how phenotypic structuring impacts the spatial and temporal dynamics of cell populations, demonstrating that different environmental pressures and phenotypic transition mechanisms significantly influence invasion patterns.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Harmful algal blooms (HABs) are characterized by rapid growth of algae, and can be caused by toxin-producing phytoplankton (TPP). When HABs occur, oxygen in the water depletes and thus can kill fish and other marine creatures causing both environmental and economic damages. In this talk, we consider a zooplankton-phytoplankton model under the assumption that the TPP exhibits group defense so that zooplankton predation decreases at high TPP density. Furthermore, we assume that toxin liberation by the TPP is not instantaneous but is rather mediated by a time lag, which is also known as the toxin liberation delay (TLD). Our results show that the model system undergoes a Hopf bifurcation around a coexistence equilibrium when the value of the TLD reaches a certain threshold. For values of the TLD just above the threshold, the stable limit cycle that is created depicts the manageable periodic fluctuation of the populations. However, when the value of the TLD is increased further, recurring blooms of various periodicity were observed which can be attributed to the occurrence of period-doubling bifurcations.

The short-tailed field vole (microtus agrestis) is the most abundant mammal in the UK, with populations reaching up to 80 million individuals. However, voles experience huge fluctuations in population numbers with up to a tenfold change over the course of regular 3.5 year cycles. Previous research has aimed to understand the mechanics behind these oscillations, but most of this work focuses on tundra regions. The ongoing Glen Finglas grazing experiment spans over 20 years, focusing on how managing grazing pressures affects various groups of species - including voles - in the more temperate upland acid grasslands of Scotland. Here, I will first present new data analysis on the Glen Finglas experiment, in particular the relationship between voles and the vegetation they use as a food source and shelter. I will then show how this data can be used to create and fit a mathematical model, capturing the vole's complex life history and interactions. Understanding these small animals is important as they are a key prey species for many predators and can also cause massive damage to plants and tree saplings. This mathematical model furthers our understanding of vole dynamics in temperate regions.

Understanding the mechanisms behind the spatial distribution, self-organisation and aggregation of organisms is a central issue in both ecology and cell biology. Since self-organisation at the population level emerges from individual behaviour, a mathematical approach is essential to elucidate these dynamics. In nature, individuals - whether cells or animals - inspect their environment before moving. This process is typically nonlocal, meaning that individuals gather information from a part of their environment rather than just their immediate location. Empirical research increasingly highlights nonlocality as a key aspect of movement, while mathematical models incorporating nonlocal interactions have gained attention for their ability to describe how interactions shape movement, reproduction and well-being. In this talk, I will present a study of a class of advection-diffusion equations that model population movement driven by nonlocal species interactions. Using a combination of analytical and numerical tools, I will show that these models support a wide range of spatio-temporal patterns, including segregation, aggregation, time-periodic behaviour, and chase-and-run dynamics. I will also discuss the existence of parameter regions with multiple stable solutions and hysteresis phenomena. Overall, I will explore various methods for analysing the bifurcation and pattern formation properties of these models, which provide essential mathematical tools for understanding the many aggregation phenomena observed in nature.

Coral reefs are some of the most threatened ecosystems on the planet but also some of the most important, hosting upwards of 25% of marine biodiversity while also providing food and livelihood to almost 1 billion people. Coral reefs are also connected together into reef networks by coral larval dispersal, meaning that management or damage at one reef may have consequences for any reefs it is connected to. For this reason, coral reef management is of interest to many industries (e.g., fisheries, tourism) and governments. The potential impact of tourism in the context of coral reef management is unclear, as tourism is a source of damage for reefs but may also be a source of income that motivates conservation actors to keep reefs healthy for future tourists. Tourism groups also often focus on a subset of coral reefs, meaning that any management initiatives driven by tourism income would also only be focused on a subset of coral reefs in a reef network. We develop a socio-ecological model composed of a system of differential equations. This model represents a network of coral reefs visited by tourists to determine whether tourism income could help sustain a healthy network of coral reefs into the future. We explore this question under a variety of different tourism paradigms, management methods (coral restoration, fisheries management) and network types. Overall, we find that management funded by tourism can help counteract tourism damage, but is unable to save reefs that are unhealthy (low initial coral cover, high fishing). Management also has limited potential to help connected reefs in the network. This study demonstrates the limited effectiveness of tourism to drive coral reef conservation and instead encourages active investment in management methods that focus on the entire reef network.

We propose a mathematical model describing the dynamics of photo fermentative bacteria leading to hydrogen production and polyhydroxybutyrate accumulation in an engineered environment. The model is derived from mass balance principles and consists of a system of differential equations describing the biomass growth, the substrate degradation and conversion into hydrogen and other catabolites, such as intracellular polyhydroxybutyrate, a precursor for bioplastics. The model accounts for crucial inhibiting phenomena and catabolic interactions affecting the evolution of the process. The study of the model has been performed also in terms of calibration with real experimental data related to specific photo-fermentative species, and it is supported by a sensitivity analysis study. The effective application of photo-fermentation for the concomitant hydrogen production and polyhydroxybutyrate accumulation was investigated.

Mosquito-borne diseases pose a significant public health challenge, and effective prevention requires accurate forecasting of mosquito populations. In this study, we developed a statistical forecasting framework that leverages climate factors, such as temperature and precipitation, to improve mosquito population predictions in Maricopa County, Arizona. Our approach combines adaptive modeling techniques and filtering methods to infer precise model parameters and address previously observed limitations, particularly the inability to capture spring dynamics. By incorporating an Ensemble Kalman Filter (EnKF) method, we estimated time-varying parameters (baseline population growth rate) and static parameters while resolving the spring problem observed in prior models. Using Generalized Additive Models (GAMs), we forecasted the baseline population growth rate on a weekly basis, integrating precipitation and temperature data as covariates. These forecasts were further used to run a mechanistic ordinary differential equation (ODE) model to predict mosquito abundance and estimate associated uncertainties. Our iterative framework was applied weekly over a 52-week period, successfully capturing seasonal variations in mosquito populations from 2014 to 2015. The EnKF demonstrated superior performance compared to traditional Markov Chain Monte Carlo (MCMC) approaches for fitting mosquito abundance data. This enhanced methodology provides actionable insights for public health decision-makers, supporting resource allocation and improving outcomes in mosquito-borne disease prevention. Our findings underscore the value of integrating climate data and adaptive filtering techniques to address forecasting challenges, ultimately enabling more effective responses to emerging or reemerging pathogens of mosquito-borne disease risks, which can be driven by human behavior to become a pandemic.

Plant-arthropod relationships are crucial to the health of global ecosystems and food production. Through co-evolution, arthropods have acquired a variety of novel senses in response to the emergence of floral cues such as scent, colour and shape. The recent discovery that several terrestrial arthropods can sense electrical fields (e-fields) motivates the investigation of floral e-fields as part of their wider sensory ecology. That is, how does a flower's morphology and material properties produce and propagate detectable, ecologically relevant electrical signals? To investigate this, we modelled the e-field interior and exterior of a flower using a novel modification of the popular AAA-least squares algorithm, extending it to two domain boundary value problems. Physically, flowers typically act as dielectrics that inductively charge in the presence of a background electrical field, e.g. charged pollinators or the Earth's atmospheric potential gradient. We therefore present the development and application of this new method for these cases and discuss the biological relevance of the results for sensory and ecological studies. Our adapted AAA algorithm gives accurate and rapid results dependent on only three parameters: the relative permittivity of the flower, flower shape and the location of the pollinator(s). The results show how flowers display distinct information about their morphology, pollen availability and nearby pollinators, at distance, through the perturbed e-field. As well as how predators, such as the crab spider, can use flowers to mask their own electrical presence and draw in unsuspecting prey. The results of the two-dimensional AAA method also shows good qualitative agreement with equivalent three-dimensional finite element models. Biologically, our results highlight the significant role floral electrics may play in plant-pollinator and predator-prey relationships, unveiling previously unstudied facets.

In this study, we utilize the species Pogonomyrmex californicus, a type of queen ant that exhibits two distinct behavioral subtypes, as inspiration for our mathematical model. The first, known as solitary queens, establish colonies independently and represent the ancestral lineage. The second subtype, cooperative queens, form groups that collectively found a single colony—an evolutionary adaptation. Laboratory experiments have revealed that these queen types display distinct behavioral traits, or personalities. To better understand the ecological implications of these differences, we develop an ordinary differential equations (ODE) model that incorporates the effects of resource availability and social organization among adult ants, particularly in relation to brood care and foraging behaviors. Our model allows for Hopf bifurcations, enabling us to analyze the conditions under which colony coexistence is promoted or collapses. With this framework, we seek to address the following key questions: How does dependency on resource availability impact the survival of each queen type? How do social organization and resource dependency together influence queen survival, and what new conditions must be met for their persistence?

Overharvesting is a pressing global problem, and spatial management, such as protecting designated areas, is one proposed solution. This talk examines how dispersal between protected and harvested areas affects the asymptotic total population size and the asymptotic yield, which are key questions for conservation management and the design of protected areas. We utilize a two-patch model with heterogeneous habitat qualities, symmetric dispersal and density-dependent growth functions in both discrete and continuous time. One patch is subject to proportional harvesting, while the other one is protected. Our results demonstrate that increased dispersal does not always increase the asymptotic total population size or the asymptotic yield. Depending on the circumstances, dispersal enables the protected patch to rescue the harvested patch from overexploitation, potentially increasing both total population size and yield. However, high levels of dispersal can also lead to a lower total population size or even cause extinction of both patches if harvesting pressure is strong. The population in the protected patch needs to have high reproductive potential and the patch needs to be the effectively larger patch in order to benefit monotonically from increased dispersal. These findings provide a fundamental understanding of how dispersal influences dynamics in fragmented landscapes under harvesting pressure.

Understanding how organisms make choices about what to eat is a fascinating puzzle explored in this study, which employs stoichiometric modeling and optimal forag- ing principles. The research delves into the intricate balance of nutrient intake with foraging strategies, investigating quality and quantity-based food selection through mathematical models. The stoichiometric models in this study, encompassing pro- ducers and a grazer, unveils the dynamics of decision-making processes, introducing fixed and variable energetic foraging costs. Analysis reveals cell quota-dependent pre- dation behaviors, elucidating biological phenomena such as “compensatory foraging behaviors” and the “stoichiometric extinction effect”. The Marginal Value Theorem quantifies food selection, highlighting the profitability of prey items and emphasizing its role in optimizing foraging strategies in predator–prey dynamics. The environ- mental factors like light and nutrient availability prove pivotal in shaping optimal foraging strategies, with numerical results from a multi-species model contributing to a comprehensive understanding of the intricate interplay between organisms and their environment.

Phototaxis, the ability of microorganisms to move in response to light, plays a crucial role in shaping the dynamics of phototrophic biofilms. While sessile cells remain typically embedded within the extracellular polymeric matrix, planktonic cells can navigate through the biofilm’s porous structure, adjusting their position in response to light cues. This directed movement optimizes exposure to favorable light conditions while avoiding harmful intensities, influencing the spatial organization and development of the biofilm community. In this talk, I will present a mathematical model for planktonic cell invasion in biofilms, where phototaxis acts as a driver of directed motility. The model incorporates a volume-filling term into the transport equation for planktonic cells, enabling the representation of phototactic behavior. A light-dependent sensitivity function captures both positive and negative phototaxis, governing cell movement toward favorable light conditions and away from excessive illumination. The biofilm is modeled as a homogeneous, viscous, incompressible fluid, with velocity described by Darcy’s law. The governing equations are solved numerically to explore the role of phototaxis in shaping biofilm dynamics. Numerical simulations reveal that motile cells accumulate in well-lit regions, enhancing sessile phototrophic growth and promoting biofilm development. The distribution of phototrophic biomass results from the interplay between random diffusion and phototactic movement. Under high-light stress conditions, photoinhibition reduces phototrophic growth and reverses phototaxis, slowing overall biofilm growth. Additionally, biofilm density modulates light penetration, either limiting phototrophic growth or providing protection against excessive exposure. These findings offer valuable insights into biofilm behavior in natural environments and can guide the optimization of biofilm-based processes in fields like wastewater treatment and bioremediation.

Evolutionary rescue occurs when a population facing environmental stress avoids extinction by rapidly acquiring beneficial mutations. While higher mutation rates can enhance rescue, the role of mutators—genotypes with elevated mutation rates—remains unclear. We develop a theoretical framework and use stochastic simulations to investigate how mutators emerge and fix under selective pressure. Our results show that mutators cannot persist in stable environments but are favored when environmental deterioration occurs, with their fixation probability influenced by the speed of environmental change and wildtype mutation rates. Pre-existing mutators further increase rescue likelihood, particularly under rapid environmental shifts. These findings provide insights into antibiotic resistance, cancer evolution, and adaptation to climate change by highlighting how environmental stress shapes mutation rate evolution.

The supplementation of additional food (AF) to introduced predators has been explored as a strategy to enhance pest control. However, AF models with prey-dependent functional responses can lead to unbounded predator growth. To address this, an AF model incorporating mutual interference has been proposed, demonstrating that pest eradication is feasible when the AF quantity ξξ exceeds a threshold function of the interference parameter ϵϵ. In this study, we revisit and extend this model, uncovering novel dynamical behaviors. We show that pest eradication occurs within a narrower AF range and can be bi-stable or globally attracting, arising through a saddle-node bifurcation. Additionally, we identify Hopf and global homoclinic bifurcations, revealing a unique dynamic where the pest extinction state becomes an 'almost' global attractor. This is the first analytical proof of such a structure in AF models, providing insights into bio-control strategies under varying predator interference conditions.

The impacts of invasive alien pests and diseases are routinely estimated and measured in the context of agriculture, but less so in the context of biodiversity and ecosystem services. In this study, we developed a new 'contribution modelling' approach to systematically estimate the impacts of pests and diseases at a continent scale. We developed this method using a case study of myrtle rust in Australia. We estimated the potential reduction of carbon sequestration in Australia due to myrtle rust using various national and scientific ecological datasets. We found that myrtle rust could lead to over 3% loss in national annual carbon sequestration if it were to spread across Australia, or over $700 million AUD value loss. While developed using a case study, this model is designed to be readily adaptible to other species and their impact on other environmental assets. Our work shows the need to systematically compile the potential impacts and costs of invasive pests and diseases to the environment and ecosystem services globally, to support both biosecurity decision-making and climate-change related initiatives such as net-zero emissions targets and reforestation efforts.

Mosquito-borne diseases pose a significant public health challenge, and effective prevention requires accurate forecasting of mosquito populations. In this study, we developed a statistical forecasting framework that leverages climate factors, such as temperature and precipitation, to improve mosquito population predictions in Maricopa County, Arizona. Our approach combines adaptive modeling techniques and filtering methods to infer precise model parameters and address previously observed limitations, particularly the inability to capture spring dynamics. By incorporating an Ensemble Kalman Filter (EnKF) method, we estimated time-varying parameters (baseline population growth rate) and static parameters while resolving the spring problem observed in prior models. Using Generalized Additive Models (GAMs), we forecasted the baseline population growth rate on a weekly basis, integrating precipitation and temperature data as covariates. These forecasts were further used to run a mechanistic ordinary differential equation (ODE) model to predict mosquito abundance and estimate associated uncertainties. Our iterative framework was applied weekly over a 52-week period, successfully capturing seasonal variations in mosquito populations from 2014 to 2015. The EnKF demonstrated superior performance compared to traditional Markov Chain Monte Carlo (MCMC) approaches for fitting mosquito abundance data. This enhanced methodology provides actionable insights for public health decision-makers, supporting resource allocation and improving outcomes in mosquito-borne disease prevention. Our findings underscore the value of integrating climate data and adaptive filtering techniques to address forecasting challenges, ultimately enabling more effective responses to emerging or reemerging pathogens of mosquito-borne disease risks, which can be driven by human behavior to become a pandemic.

Plant-arthropod relationships are crucial to the health of global ecosystems and food production. Through co-evolution, arthropods have acquired a variety of novel senses in response to the emergence of floral cues such as scent, colour and shape. The recent discovery that several terrestrial arthropods can sense electrical fields (e-fields) motivates the investigation of floral e-fields as part of their wider sensory ecology. That is, how does a flower's morphology and material properties produce and propagate detectable, ecologically relevant electrical signals? To investigate this, we modelled the e-field interior and exterior of a flower using a novel modification of the popular AAA-least squares algorithm, extending it to two domain boundary value problems. Physically, flowers typically act as dielectrics that inductively charge in the presence of a background electrical field, e.g. charged pollinators or the Earth's atmospheric potential gradient. We therefore present the development and application of this new method for these cases and discuss the biological relevance of the results for sensory and ecological studies. Our adapted AAA algorithm gives accurate and rapid results dependent on only three parameters: the relative permittivity of the flower, flower shape and the location of the pollinator(s). The results show how flowers display distinct information about their morphology, pollen availability and nearby pollinators, at distance, through the perturbed e-field. As well as how predators, such as the crab spider, can use flowers to mask their own electrical presence and draw in unsuspecting prey. The results of the two-dimensional AAA method also shows good qualitative agreement with equivalent three-dimensional finite element models. Biologically, our results highlight the significant role floral electrics may play in plant-pollinator and predator-prey relationships, unveiling previously unstudied facets.

In this study, we utilize the species Pogonomyrmex californicus, a type of queen ant that exhibits two distinct behavioral subtypes, as inspiration for our mathematical model. The first, known as solitary queens, establish colonies independently and represent the ancestral lineage. The second subtype, cooperative queens, form groups that collectively found a single colony—an evolutionary adaptation. Laboratory experiments have revealed that these queen types display distinct behavioral traits, or personalities. To better understand the ecological implications of these differences, we develop an ordinary differential equations (ODE) model that incorporates the effects of resource availability and social organization among adult ants, particularly in relation to brood care and foraging behaviors. Our model allows for Hopf bifurcations, enabling us to analyze the conditions under which colony coexistence is promoted or collapses. With this framework, we seek to address the following key questions: How does dependency on resource availability impact the survival of each queen type? How do social organization and resource dependency together influence queen survival, and what new conditions must be met for their persistence?

Overharvesting is a pressing global problem, and spatial management, such as protecting designated areas, is one proposed solution. This talk examines how dispersal between protected and harvested areas affects the asymptotic total population size and the asymptotic yield, which are key questions for conservation management and the design of protected areas. We utilize a two-patch model with heterogeneous habitat qualities, symmetric dispersal and density-dependent growth functions in both discrete and continuous time. One patch is subject to proportional harvesting, while the other one is protected. Our results demonstrate that increased dispersal does not always increase the asymptotic total population size or the asymptotic yield. Depending on the circumstances, dispersal enables the protected patch to rescue the harvested patch from overexploitation, potentially increasing both total population size and yield. However, high levels of dispersal can also lead to a lower total population size or even cause extinction of both patches if harvesting pressure is strong. The population in the protected patch needs to have high reproductive potential and the patch needs to be the effectively larger patch in order to benefit monotonically from increased dispersal. These findings provide a fundamental understanding of how dispersal influences dynamics in fragmented landscapes under harvesting pressure.

Understanding how organisms make choices about what to eat is a fascinating puzzle explored in this study, which employs stoichiometric modeling and optimal forag- ing principles. The research delves into the intricate balance of nutrient intake with foraging strategies, investigating quality and quantity-based food selection through mathematical models. The stoichiometric models in this study, encompassing pro- ducers and a grazer, unveils the dynamics of decision-making processes, introducing fixed and variable energetic foraging costs. Analysis reveals cell quota-dependent pre- dation behaviors, elucidating biological phenomena such as “compensatory foraging behaviors” and the “stoichiometric extinction effect”. The Marginal Value Theorem quantifies food selection, highlighting the profitability of prey items and emphasizing its role in optimizing foraging strategies in predator–prey dynamics. The environ- mental factors like light and nutrient availability prove pivotal in shaping optimal foraging strategies, with numerical results from a multi-species model contributing to a comprehensive understanding of the intricate interplay between organisms and their environment.

Phototaxis, the ability of microorganisms to move in response to light, plays a crucial role in shaping the dynamics of phototrophic biofilms. While sessile cells remain typically embedded within the extracellular polymeric matrix, planktonic cells can navigate through the biofilm’s porous structure, adjusting their position in response to light cues. This directed movement optimizes exposure to favorable light conditions while avoiding harmful intensities, influencing the spatial organization and development of the biofilm community. In this talk, I will present a mathematical model for planktonic cell invasion in biofilms, where phototaxis acts as a driver of directed motility. The model incorporates a volume-filling term into the transport equation for planktonic cells, enabling the representation of phototactic behavior. A light-dependent sensitivity function captures both positive and negative phototaxis, governing cell movement toward favorable light conditions and away from excessive illumination. The biofilm is modeled as a homogeneous, viscous, incompressible fluid, with velocity described by Darcy’s law. The governing equations are solved numerically to explore the role of phototaxis in shaping biofilm dynamics. Numerical simulations reveal that motile cells accumulate in well-lit regions, enhancing sessile phototrophic growth and promoting biofilm development. The distribution of phototrophic biomass results from the interplay between random diffusion and phototactic movement. Under high-light stress conditions, photoinhibition reduces phototrophic growth and reverses phototaxis, slowing overall biofilm growth. Additionally, biofilm density modulates light penetration, either limiting phototrophic growth or providing protection against excessive exposure. These findings offer valuable insights into biofilm behavior in natural environments and can guide the optimization of biofilm-based processes in fields like wastewater treatment and bioremediation.

Evolutionary rescue occurs when a population facing environmental stress avoids extinction by rapidly acquiring beneficial mutations. While higher mutation rates can enhance rescue, the role of mutators—genotypes with elevated mutation rates—remains unclear. We develop a theoretical framework and use stochastic simulations to investigate how mutators emerge and fix under selective pressure. Our results show that mutators cannot persist in stable environments but are favored when environmental deterioration occurs, with their fixation probability influenced by the speed of environmental change and wildtype mutation rates. Pre-existing mutators further increase rescue likelihood, particularly under rapid environmental shifts. These findings provide insights into antibiotic resistance, cancer evolution, and adaptation to climate change by highlighting how environmental stress shapes mutation rate evolution.

The supplementation of additional food (AF) to introduced predators has been explored as a strategy to enhance pest control. However, AF models with prey-dependent functional responses can lead to unbounded predator growth. To address this, an AF model incorporating mutual interference has been proposed, demonstrating that pest eradication is feasible when the AF quantity ξξ exceeds a threshold function of the interference parameter ϵϵ. In this study, we revisit and extend this model, uncovering novel dynamical behaviors. We show that pest eradication occurs within a narrower AF range and can be bi-stable or globally attracting, arising through a saddle-node bifurcation. Additionally, we identify Hopf and global homoclinic bifurcations, revealing a unique dynamic where the pest extinction state becomes an 'almost' global attractor. This is the first analytical proof of such a structure in AF models, providing insights into bio-control strategies under varying predator interference conditions.

The impacts of invasive alien pests and diseases are routinely estimated and measured in the context of agriculture, but less so in the context of biodiversity and ecosystem services. In this study, we developed a new 'contribution modelling' approach to systematically estimate the impacts of pests and diseases at a continent scale. We developed this method using a case study of myrtle rust in Australia. We estimated the potential reduction of carbon sequestration in Australia due to myrtle rust using various national and scientific ecological datasets. We found that myrtle rust could lead to over 3% loss in national annual carbon sequestration if it were to spread across Australia, or over $700 million AUD value loss. While developed using a case study, this model is designed to be readily adaptible to other species and their impact on other environmental assets. Our work shows the need to systematically compile the potential impacts and costs of invasive pests and diseases to the environment and ecosystem services globally, to support both biosecurity decision-making and climate-change related initiatives such as net-zero emissions targets and reforestation efforts.

Climate change is altering the mean, variance and autocorrelation of temperature. However, linear approaches to incorporating these temperature impacts in simple population models do not provide realistic predictions regarding climate change impacts. For example, simple degree days approaches or using a linear function of temperature to alter the density-independent population growth rate will not account for the sometimes catastrophic decrease in performance with high temperatures. We use an extremely simple population model coupled to nonlinear thermal performance curves to explore the simultaneous impact of changes to temperature mean, variance and autocorrelation. The realized density-independent population growth rate is given by three types of thermal performance curves that correspond to published data. We find relatively small impacts on established population dynamics when realistic changes in temperature sequences are used, suggesting that many populations may be quite robust to temperature-driven climate change impacts in the near term. The most extreme right-skewed performance curves are most likely to result in species extinctions, even though these curves have higher optimal temperatures.

The oriental fruit fly, Bactrocera dorsalis, is a serious threat to crops and orchards in many places around the World, and in particular in Réunion island, where it was first detected in 2017. Since then, this pest has invaded the whole island and displaced established fruit fly populations. Since Réunion island is a hot spot of diversity, appropriate control tools have to be deployed to eliminate or reduce the wild population. I will present recent results that study the combination of the Sterile Insect Technique, entomopathogen fungi, and also pheromone traps. In particular, we will show how the spatial component and the orchards connectivity can drastically change the releases strategy, as well as the critical amount of sterile insects to release. We discuss (optimal) strategies obtained with our models versus realistic strategies that can actually be developed in the field. Our approach being generic, it can be adapted to other pests and disease vectors, such as mosquitoes. This works stands within the AttracTIS project, funded by Ecophyto 2021-2022.

As parasitoids are the most frequently used biocontrol agents, especially in agriculture and forest ecosystems, they have become a cornerstone in mathematical biology. They are also a prototypical example of discrete-time systems. Here we consider a simple host-parasitoid model that is based on the classical Nicholson-Bailey model, but includes two extensions that are ecologically plausible: (1) density-dependent host growth (of Beverton-Holt type) and (2) a functional response of type III. The latter can be caused by a number of ecological mechanisms and is key in driving a rich dynamical behavior. While the system admits at most one nontrivial fixed point, we observe up to four coexisting non-equilibrium attractors. They can be periodic, quasi-periodic, or chaotic. They emerge in a quasi-periodic route to chaos and exhibit frequency-locking phenomena. We find different regular organized structures in the two-dimensional parameter plane that describe periodic oscillations surrounded by chaos. Among these structures are Arnold tongues (which have been previously reported in related models) and shrimp-shaped domains, which are little known in ecological models. Our results demonstrate that a type III functional response of parasitoids induces many new complex phenomena. While in continuous-time models the type III functional response tends to be stabilizing, in discrete-time models it can have very contrasting effects. The ecological implications are a high sensitivity not only to parameters but also to the initial condition.

The site frequency spectrum (SFS) is a popular summary statistic of genomic data. In population genetics, the SFS has provided a simple means of inferring the rate of adaptation of a population and for distinguishing between neutral evolution and evolution under selection. The rapidly growing amount of cancer genomic data has attracted interest in the SFS of an exponentially growing population. In this talk, we discuss recent results on the expected value of the SFS of a population that grows according to a stochastic branching process, as well as (first-order) almost sure convergence results for the SFS in the large-time and large-population limits. Our results show that while the SFS depends linearly on the mutation rate, the branching process parameters of birth and death control the fundamental shape of the SFS at the low-frequency end. For the special case of a birth-death process (binary branching process), our results give rise to statistically consistent estimators for the mutation rate and extinction probability of the population, which stands in contrast to previous work which has indicated the need for additional data to decouple these two parameters. Overall, our work shows how single timepoint data on the SFS of an exponentially growing population can be used to infer important evolutionary parameters.

Wildfires disrupt ecosystems, with climate change exacerbating vulnerability in regions poorly adapted to such disturbances. These events are driven by complex, multi-scale interactions where small perturbations in environmental factors can trigger large-scale shifts, complicating prediction efforts. We propose a coupled convection-reaction-diffusion system as a framework for modelling wildfire spread dynamics. This system integrates spatial and temporal variability to identify thresholds for spread and quantify the impact of abrupt environmental changes on burnt areas and rates of propagation. Incorporating environmental, meteorological, and historical fire record data from the Global Wildfire Information System, the Department for Environment, Food and Rural Affairs (UK), and drone footage of heather burning. Bayesian inference and Monte Carlo methods are employed for parameter estimation and uncertainty quantification, ensuring robust model validation against unseen data. Recent wildfire events around the globe highlight the need for actionable insights into environmental vulnerability, property loss, and infrastructure risk. By enabling near-real-time simulations, this model aims to provide a computational tool for emergency response, long-term management strategies, and assessments of climate change-induced outlier weather patterns influencing fire behaviour. This work highlights the potential of mathematical modelling to advance understanding and management of critical ecological disturbances.

Two fundamental quantities in population biology, the reproductive number R0 and the growth rate, are intimately linked, but the exact nature of their relationship is somewhat obscure. Models of microbial growth typically have R0=2, but estimating their growth rate, and hence fitness, requires solving the famous Euler-Lotka equation. Conversely, in epidemiology one typically measures how quickly the infected population grows, but it is the reproductive number R0 that sets the threshold for an epidemic breakout and for herd immunity. In this talk, we use statistical techniques based on large deviations theory to clarify how exactly the population growth rate and R0 are connected. Building an analogy to classical thermodynamics, we show that the long-term behaviour of a population is encoded in a single convex function that relates growth rate, R0, and the statistics of intergeneration times in lineages. As an application, we derive a general formulation of the Euler-Lotka equation and explain why it is almost always appears as an implicit equation.

In the 90s, several experiments suggested a hypothesis that certain genes function to mask or buffer the effects of mutations, thereby allowing them to accumulate and be stored. These were termed evolutionary capacitors and addressed the fundamental evolutionary problem of how populations optimize fitness in one environment while maintaining variation to adapt to another. However, more recent experiments support an alternative hypothesis that such buffering of mutations is a natural and unsurprising outcome of epistasis and the mutation-selection process. To quantitatively test this hypothesis, we develop a mathematical framework that extends a classical partial differential equation of the mutation-selection process to account for epistasis. Using a perturbation method on steady state solutions, we show that certain types of epistatic interactions and selection pressures will lead to the emergence of the evolutionary capacitance phenomena.

Methane, a potent greenhouse gas (GHG), is now driving climate change at an unprecedented rate. With a warming potential greater than carbon dioxide, it poses a substantial threat to the functioning of ecosystems. Despite its importance, studies investigating its direct impact on species interactions within ecosystems are rare. This growing concern highlights the need for a comprehensive understanding of the factors that could disrupt food chains, ultimately impacting ecosystem stability and resilience. In this talk, I will address this gap by developing a mechanistic model that integrates methane dynamics with the populations of species and detritus. This novel approach offers a framework for understanding how gaseous pollutants like methane influence trophic interactions. The model is studied for a range of concentrations of methane. Our findings reveal that low concentrations of methane can benefit species growth as an alternative carbon source. However, moderate to high levels induce sub-lethal to lethal effects. Further, analyzing the mechanisms for long transients in the fast-intermediate-slow formulation of the model, I will discuss how faster methane accumulation in water can result in slower species growth.

Microbially induced corrosion (MIC) is a significant global issue impacting infrastructure, economies, and environment. In wastewater systems, MIC is primarily associated with biofilm formation on concrete sewer pipes, leading to severe degradation due to microbial metabolic activity. The proliferation of sewer biofilms occurs in both submerged and unsubmerged conditions, leading to distinct microbial communities. Commonly, these biofilms host microorganisms such as fermentation bacteria, hydrogen-producing acetogens, denitrifying bacteria, sulfate-reducing bacteria, sulfur-oxidizing bacteria, and methanogens. In particular, sulfur-oxidizing bacteria play a crucial role in corrosion, as they oxidize hydrogen sulfide from wastewater effluents, generating sulfuric acid that accelerates concrete deterioration. A one-dimensional model with double free boundaries has been developed to investigate the proliferation of biofilms and the related corrosion process in wastewater concrete pipes. The domain is composed of two free boundary regions: a biofilm that grows towards the interior cavity of the pipe, sitting on a gypsum layer formed by corrosion, which penetrates the concrete pipe. Diffusion-reaction equations govern the transport and the metabolic production or consumption of dissolved substances, such as hydrogen sulfide, oxygen, and sulfuric acid within the biofilm layer. The biofilm free boundary tracks the growth of the microbial community, regulated by microbial metabolic activity and detachment phenomena. The corrosion process is incorporated into the model through a Stefan-type condition, which drives the advancement of the gypsum free boundary into the concrete pipe, governed by microbial production of sulfuric acid. Numerical simulations have been carried out to investigate the model behavior, encompassing the development and progression of the biofilm as well as the corrosion advancement, with the aim of elucidating the key factors governing both phenomena.

The present paper investigates the dynamical behavior of the modified May Holling Tanner model in the presence of dynamic alternative resources. We study the role of dynamic alternative resources on the survival of the species when there is prey rarity. Detailed mathematical analysis and numerical evaluations, including the situation of ecosystem collapsing, have been presented to discuss the coexistence of species', stability, occurrence of different bifurcations (saddle-node, transcritical, and Hopf) in three cases in the presence of prey and alternative resources, in the absence of prey and in the absence of alternative resources. It has been obtained that the multiple coexisting states and their stability are outcomes of variations in predation rate for alternative resources. Also, the occurrence of Hopf bifurcation, saddle-node bifurcation, and transcritical bifurcation are due to variations in the parameters of dynamic alternative resources. The impact of dynamic alternative resources on species' density reveals the fact that if the predation rate for alternative resources increases, then the prey biomass increases (under some restrictions), and variations in the predator's biomass widely depend upon the quality of food items. This study also points out that the survival of predators is possible in the absence of prey. In the theme of ecological balance, the present study suggests some theoretical points of view for the eco-managers.

The decision of an individual or group to leave its current environment may be influenced by various factors. These include external or inter-species factors such as the presence of predators or food availability, and also intra-species dynamics like mate searching or the strength of social ties within a group. Understanding the consequences of these behaviours on the population level dynamics is non-trivial. In this study, we explore a stochastic model describing the movement of males and females of a species between localised patches, in which the movement rates are dependent on the sex ratio within the patch. By deriving a system of stochastic differential equations governing the fluctuations in these patches we can model a diverse range of intra-species behaviours driven solely by an individual's response to local sex ratio. We subsequently uncover and explore how different individual behaviours can give rise to large scale (anti-)correlated movements between the sexes.

To comprehend the mechanisms driving biodiversity and ecosystem resilience in a rapidly changing world, it is essential to explore the behavioral diversity among individuals in greater depth. Consistent individual differences in behavior, often referred to as animal personality, play a crucial role in shaping ecological and evolutionary dynamics, particularly in foraging behavior. Traditional approaches in behavioral and evolutionary ecology typically focus on average behavior, neglecting the significance of individual variability. This study explores the influence of consumer personality on ecological dynamics, specifically examining how variations in food availability affect behavioral strategies and ecosystem functioning. We develop a resource-consumer model that incorporates personality-dependent saturating attack rates based on the mean-field ratio of resources to consumers. The well-posedness of the model is established, and we analyze the existence and stability of all steady-state solutions. Through bifurcation analysis, we identify critical transition parameters and describe the nonlinear phenomena induced by personality-dependent attack rates. Our findings demonstrate that boldness in consumers enhances their persistence, particularly under low levels of boldness, where populations can survive even with moderate or high food supply, which was not captured in classical frameworks.

A microscale ecological study using the closure approach to understand the impact of productivity controlled by geographical and seasonal variations on the intermittency of phytoplankton will be presented in this talk. Using this approach for a nutrient–phytoplankton model with Holling type III functional response, it has been shown how the dynamics of the system can be affected by the environmental fluctuations triggered by the impact of light, temperature, and salinity, which fluctuate with regional and seasonal variations. Reynold’s averaging method in space, which results in expressing the original components in terms of its mean (average value) and perturbation (fluctuation) has been used to determine the impact of growth fluctuation in phytoplankton distribution and in the intermittency of phytoplankton spreading (variance). Parameters are estimated from the nature of productivity and spread of phytoplankton density during field observation done at four different locations of Tokyo Bay. The model validation shows that our results are in good agreement with the field observation and succeeded in explaining the intermittent phytoplankton distribution at different locations of Tokyo Bay, Japan, and its neighboring coastal regions.

Two fundamental quantities in population biology, the reproductive number R0 and the growth rate, are intimately linked, but the exact nature of their relationship is somewhat obscure. Models of microbial growth typically have R0=2, but estimating their growth rate, and hence fitness, requires solving the famous Euler-Lotka equation. Conversely, in epidemiology one typically measures how quickly the infected population grows, but it is the reproductive number R0 that sets the threshold for an epidemic breakout and for herd immunity. In this talk, we use statistical techniques based on large deviations theory to clarify how exactly the population growth rate and R0 are connected. Building an analogy to classical thermodynamics, we show that the long-term behaviour of a population is encoded in a single convex function that relates growth rate, R0, and the statistics of intergeneration times in lineages. As an application, we derive a general formulation of the Euler-Lotka equation and explain why it is almost always appears as an implicit equation.

In the 90s, several experiments suggested a hypothesis that certain genes function to mask or buffer the effects of mutations, thereby allowing them to accumulate and be stored. These were termed evolutionary capacitors and addressed the fundamental evolutionary problem of how populations optimize fitness in one environment while maintaining variation to adapt to another. However, more recent experiments support an alternative hypothesis that such buffering of mutations is a natural and unsurprising outcome of epistasis and the mutation-selection process. To quantitatively test this hypothesis, we develop a mathematical framework that extends a classical partial differential equation of the mutation-selection process to account for epistasis. Using a perturbation method on steady state solutions, we show that certain types of epistatic interactions and selection pressures will lead to the emergence of the evolutionary capacitance phenomena.

Methane, a potent greenhouse gas (GHG), is now driving climate change at an unprecedented rate. With a warming potential greater than carbon dioxide, it poses a substantial threat to the functioning of ecosystems. Despite its importance, studies investigating its direct impact on species interactions within ecosystems are rare. This growing concern highlights the need for a comprehensive understanding of the factors that could disrupt food chains, ultimately impacting ecosystem stability and resilience. In this talk, I will address this gap by developing a mechanistic model that integrates methane dynamics with the populations of species and detritus. This novel approach offers a framework for understanding how gaseous pollutants like methane influence trophic interactions. The model is studied for a range of concentrations of methane. Our findings reveal that low concentrations of methane can benefit species growth as an alternative carbon source. However, moderate to high levels induce sub-lethal to lethal effects. Further, analyzing the mechanisms for long transients in the fast-intermediate-slow formulation of the model, I will discuss how faster methane accumulation in water can result in slower species growth.

Microbially induced corrosion (MIC) is a significant global issue impacting infrastructure, economies, and environment. In wastewater systems, MIC is primarily associated with biofilm formation on concrete sewer pipes, leading to severe degradation due to microbial metabolic activity. The proliferation of sewer biofilms occurs in both submerged and unsubmerged conditions, leading to distinct microbial communities. Commonly, these biofilms host microorganisms such as fermentation bacteria, hydrogen-producing acetogens, denitrifying bacteria, sulfate-reducing bacteria, sulfur-oxidizing bacteria, and methanogens. In particular, sulfur-oxidizing bacteria play a crucial role in corrosion, as they oxidize hydrogen sulfide from wastewater effluents, generating sulfuric acid that accelerates concrete deterioration. A one-dimensional model with double free boundaries has been developed to investigate the proliferation of biofilms and the related corrosion process in wastewater concrete pipes. The domain is composed of two free boundary regions: a biofilm that grows towards the interior cavity of the pipe, sitting on a gypsum layer formed by corrosion, which penetrates the concrete pipe. Diffusion-reaction equations govern the transport and the metabolic production or consumption of dissolved substances, such as hydrogen sulfide, oxygen, and sulfuric acid within the biofilm layer. The biofilm free boundary tracks the growth of the microbial community, regulated by microbial metabolic activity and detachment phenomena. The corrosion process is incorporated into the model through a Stefan-type condition, which drives the advancement of the gypsum free boundary into the concrete pipe, governed by microbial production of sulfuric acid. Numerical simulations have been carried out to investigate the model behavior, encompassing the development and progression of the biofilm as well as the corrosion advancement, with the aim of elucidating the key factors governing both phenomena.

The present paper investigates the dynamical behavior of the modified May Holling Tanner model in the presence of dynamic alternative resources. We study the role of dynamic alternative resources on the survival of the species when there is prey rarity. Detailed mathematical analysis and numerical evaluations, including the situation of ecosystem collapsing, have been presented to discuss the coexistence of species', stability, occurrence of different bifurcations (saddle-node, transcritical, and Hopf) in three cases in the presence of prey and alternative resources, in the absence of prey and in the absence of alternative resources. It has been obtained that the multiple coexisting states and their stability are outcomes of variations in predation rate for alternative resources. Also, the occurrence of Hopf bifurcation, saddle-node bifurcation, and transcritical bifurcation are due to variations in the parameters of dynamic alternative resources. The impact of dynamic alternative resources on species' density reveals the fact that if the predation rate for alternative resources increases, then the prey biomass increases (under some restrictions), and variations in the predator's biomass widely depend upon the quality of food items. This study also points out that the survival of predators is possible in the absence of prey. In the theme of ecological balance, the present study suggests some theoretical points of view for the eco-managers.

The decision of an individual or group to leave its current environment may be influenced by various factors. These include external or inter-species factors such as the presence of predators or food availability, and also intra-species dynamics like mate searching or the strength of social ties within a group. Understanding the consequences of these behaviours on the population level dynamics is non-trivial. In this study, we explore a stochastic model describing the movement of males and females of a species between localised patches, in which the movement rates are dependent on the sex ratio within the patch. By deriving a system of stochastic differential equations governing the fluctuations in these patches we can model a diverse range of intra-species behaviours driven solely by an individual's response to local sex ratio. We subsequently uncover and explore how different individual behaviours can give rise to large scale (anti-)correlated movements between the sexes.

To comprehend the mechanisms driving biodiversity and ecosystem resilience in a rapidly changing world, it is essential to explore the behavioral diversity among individuals in greater depth. Consistent individual differences in behavior, often referred to as animal personality, play a crucial role in shaping ecological and evolutionary dynamics, particularly in foraging behavior. Traditional approaches in behavioral and evolutionary ecology typically focus on average behavior, neglecting the significance of individual variability. This study explores the influence of consumer personality on ecological dynamics, specifically examining how variations in food availability affect behavioral strategies and ecosystem functioning. We develop a resource-consumer model that incorporates personality-dependent saturating attack rates based on the mean-field ratio of resources to consumers. The well-posedness of the model is established, and we analyze the existence and stability of all steady-state solutions. Through bifurcation analysis, we identify critical transition parameters and describe the nonlinear phenomena induced by personality-dependent attack rates. Our findings demonstrate that boldness in consumers enhances their persistence, particularly under low levels of boldness, where populations can survive even with moderate or high food supply, which was not captured in classical frameworks.

A microscale ecological study using the closure approach to understand the impact of productivity controlled by geographical and seasonal variations on the intermittency of phytoplankton will be presented in this talk. Using this approach for a nutrient–phytoplankton model with Holling type III functional response, it has been shown how the dynamics of the system can be affected by the environmental fluctuations triggered by the impact of light, temperature, and salinity, which fluctuate with regional and seasonal variations. Reynold’s averaging method in space, which results in expressing the original components in terms of its mean (average value) and perturbation (fluctuation) has been used to determine the impact of growth fluctuation in phytoplankton distribution and in the intermittency of phytoplankton spreading (variance). Parameters are estimated from the nature of productivity and spread of phytoplankton density during field observation done at four different locations of Tokyo Bay. The model validation shows that our results are in good agreement with the field observation and succeeded in explaining the intermittent phytoplankton distribution at different locations of Tokyo Bay, Japan, and its neighboring coastal regions.

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

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

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

What should you do if you miss a dose of antibiotics? Despite the prevalence of missed antibiotic doses, there is vague or little guidance on what to do when a dose is forgotten. In this paper, we consider the effects of different patient responses after missing a dose using a mathematical model that links antibiotic concentration with bacteria dynamics. We show using simulations that, in some circumstances, (a) missing just a few doses can cause treatment failure, and (b) this failure can be remedied by simply taking a double dose after a missed dose. We then develop an approximate model that is analytically tractable and use it to understand when it might be advisable to take a double dose after a missed dose.

Natural killer (NK) cells are members of the innate immune system and are proving to be a lethal weapon against cancer. To unlock the full power of NK cells, we must first address the central question: How does an NK cell recognize a malignant cell? To assess a target cell, an NK cell forms an immunological synapse, which is the interaction zone between the two cells. Ligand-receptor binding within the synapse triggers downstream activating and inhibitory signaling pathways that integrate to control the actin cytoskeleton network. Dominating activating signals causes the NK cell’s actin network to reorganize which transports more receptors to the synapse, thereby generating a positive feedback loop. Mechanistically, activating signals lead to the activation of the Arp2/3 complex which creates a branched actin network. In return, the flow of this network drives the centripetal transport of receptors to the synapse. We propose an advection-diffusion model to capture this phenomenon. Furthermore, we test what ligand-receptor densities permit synapse formation.

Kaposi’s sarcoma-associated herpesvirus (KSHV) is a causative agent of several lymphoproliferative diseases, particularly in immunocompromised individuals. These malignancies originate from latently infected B cells, where KSHV persists as extrachromosomal episomes. While the viral protein LANA is essential for viral maintenance during latency, the mechanisms enabling lifelong persistence remain unclear. To quantify episome dynamics, we developed a mathematical model of latent KSHV replication and segregation during cell division, and a statistical framework to infer viral dynamics from fluorescent microscopy images. We built a Gibbs sampler to extract episome counts from imperfectly resolved images of pre- and post-division cells. Using these counts, we estimate the efficiency of replication and segregation, propagating imaging uncertainty into our parameter estimates. Our framework, validated on synthetic data, provided similar estimates of replication efficiency (78%, 95% CI [53%, 90%]) and segregation efficiency (91% [78%, 100%]) when applied to fixed and live images of cells transfected with either full-length KSHV or a minimal plasmid capable of episome maintenance. Simulations of a dividing cell population showed that imperfect replication and segregation preclude decades-long persistence without the assistance of additional mechanisms such as cell-survival benefits to infection or occasional lytic replication. We also modeled KSHV-dependent malignancies to evaluate episome replication and segregation as targets to control tumor growth. Simulations revealed that reducing replication effectively disrupts tumor growth, with the required reduction dependent on cell division kinetics. Our results suggest that KSHV employs a partitioning mechanism, as opposed to random segregation, though replication and segregation are imperfect. Furthermore, targeting episome replication may offer a viable strategy to reduce tumor burden in KSHV-associated malignancies.

Individual bacterium make decisions at a genetic level as a result of various types of gene regulation; this process plays out on a population level to inform colony growth. Vibrio vulnificus is an opportunistic Gramnegative marine pathogen with a limiting growth factor of iron. Compared to other foodborne pathogens, Vibrio vulnificus has a high mortality rate and relatively poorly understood virulence mechanisms. When inside a human host, this bacteria utilizes heme as a source of iron, necessitating the ability to turn pieces of the heme acquisition system off and on in response to various environmental signals. As establishment of infection depends on Vibrio vulnificus’s ability to change from a marine to human environment, the ability to switch on the heme-intake system is an important part of establishment of initial infection. One such part of this system is the hemolysin VvhA. This toxin is excreted by the bacterium to lyse erythrocytes, thereby releasing heme into the extracellular environment where the bacteria can use it as a source of iron. This toxin is regulated by a complex set of factors including nutrient availability and quorum sensing. Exploring this gene regulatory network via bifurcation analysis reveals a complex bifurcation structure. These dynamics allow an individual bacterium to integrate a variety of signals in response to a changing environment. In particular, bistability in the system points to the likelihood of a heterogenous bacterial colony, where many bacteria benefit from a smaller number of hemolysin producers. This allows for modeling both a heterogeneous population and incorporation of the physiological mechanism by which cells make the decision to switch states. The interdependence between toxin production, nutrient availability, and colony growth result in interplay between the bacteria and their environment, allowing for insights into the overall course of infection.

Pertussis (whooping cough), caused by Bordetella pertussis, is most severe in infants, with most deaths occurring in unvaccinated infants under three months of age. Vaccination with the DTaP (priming) and TdaP (booster) immunizations is effective, with TdaP given during pregnancy and DTaP in infancy. However, immunomodulation can dampen the IgG response in infants born to vaccinated mothers. We hypothesize that adjusting the vaccination schedule could reduce immunomodulation and enhance vaccine efficacy. Since empirically testing multiple schedules is impractical, we propose using mathematical modeling alongside two experimental mouse models to determine an optimal schedule. Pregnant and infant mice will be immunized following a murine analog of standard vaccination. These data will inform our model, allowing us to estimate key immune parameters. Once parametrized, our model will propose schedules that maximize infant antibody response. A second mouse experiment will test these schedules, comparing immune responses to assess their efficacy. This approach will help evaluate immunomodulation mechanisms and refine vaccination strategies. The mechanistic evaluation of immunomodulation is of significance due to its lack of effective investigations to date.

In this talk, we will present effective antibiotic regimens in the presence of drug-resistant bacteria in biofilms. We begin by discussing models of in-vivo antimicrobial resistance transfer within bacterial biofilms, focusing on various one-dimensional biofilm models. Our approach includes modeling resistance acquisition through horizontal gene transfer between resistant and susceptible strains while also accounting for the role of persistor cells. We examine the effects of periodic antibiotic dosing at a constant level, showing that it may not always lead to bacterial eradication. To address this challenge, we utilize a numerical optimization algorithm to determine the optimal antibiotic dosing strategy. Additionally, we analyze how changes in different model parameters impact the qualitative behavior of the optimal dosing regimen.

Axial spondyloarthritis (SpA) is a chronic inflammatory disease impacting the joints in the axial skeleton (e.g. chest, spine, pelvis). Tumor necrosis factor inhibitors (TNFi), such as the monoclonal antibody adalimumab, are used to treat severe cases, but therapy is often discontinued due to loss of efficacy, not least resulting from immunogeniocity and development of anti-drug antibodies (ADA). The disease-modifying drug methotrexate (MTX) has shown potential in reducing the formation of ADA to various TNFi in rheumatoid arthritis (Krickaert et al, 2012), but little is known about its mode of action. We adapt a mechanistic mathematical model for immunogencity towards adalimumab based on Chen et al. 2014 to describe the impact of MTX in reducing immunogenicity, and parametrise it based on patient data from a multicentric randomised trial (Ducoureau et al, 2020). ODEs describe the pharmacokinetics and pharmacodynamics of the therapeutic compounds (adalimumab only or adalimumab and MTX, depending on patient cohort), the dynamics of T and B lymphocytes, antigen presenting cells, and some relevant cytokines for the disease. Due to the large size of this model, we employ several reduced models to estimate some of the parameter values. The model is used to simulate several scenarios in order to elucidate the most likely modes of action of MTX to reduce immunogenicity by comparing the simulated and measured ADA titres along 5 hospital visits. This is joint work with Sara Sottile (Bologna, Italy) and Denis Mulleman (Tours, France). This work is based upon work from COST Action ENOTTA (CA21147), supported by COST (European Cooperation in Science and Technology), and Contract KP-06-DKOST-13 of the Bulgarian Fund for Scientific Research.

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

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

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

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

Identifying early warning signals (EWS) of shifts in vaccinating behaviors can be helpful in predicting disease outbreaks. Evolutionary game theory has been used to model individual vaccination decisions, while bifurcation theory has identified statistical EWS, such as increasing variance and lag-1 autocorrelation, near critical transitions. However, these conventional methods often struggle with noisy data. In this study, we improve coupled behavior-disease models by incorporating population heterogeneity, distinguishing between social media users and non-users, and examining the role of homophily in shaping disease dynamics. We develop deep learning classifiers, including Long Short-Term Memory (LSTM) and Residual Neural Networks (ResNet), trained on simulated data from a stochastic coupled model with Lévy noise that captures the heavy-tailed fluctuations characteristic of real-world systems. Our results show that these models outperform traditional statistical indicators in both sensitivity and specificity while offering clearer interpretability on empirical data. This approach provides a robust framework for detecting EWS and improving outbreak prediction, highlighting the power of deep learning in real-time public health monitoring.

Evaluating vaccine efficacy (VE) during emerging epidemics is challenging due to unpredictable transmission dynamics. An age-structured SEIAR compartmental model was developed using South Korea’s 2022 population and parameters from COVID-19 and the 2009 H1N1 pandemic to optimize RCT timing and sample size. Simulations varied trial initiation (±10%, ±20%, ±30% of the epidemic peak), follow-up (4–12 weeks), recruitment (2–12 weeks), and non-pharmaceutical interventions (10–20%). Results showed that VE remained stable, but sample size requirements fluctuated, decreasing post-peak before rising sharply. Starting trials 30% before the peak with extended recruitment minimized sample sizes without compromising power. NPIs expanded trial feasibility, and sample size estimates from simulated placebo cases maintained >85% power, avoiding under- or over-powering. This model provides a framework for designing adaptive and efficient vaccine trials in future respiratory epidemics.

Social contact patterns are crucial for understanding infectious disease transmission, but detailed data has been scarce in South Korea. We conducted a two-week survey covering various periods, including school vacations and holidays. Participants provided information on their contacts, including location, duration, frequency, and characteristics of the contact person. Analysis of the data from 1,987 participants revealed 133,776 contacts, averaging 4.81 contacts per person daily. Contact numbers varied by age, household size, and time period. The highest number of contacts was observed in the 5-19 age group, lowest in the 20-29 group, and gradually increased up to the 70+ group. Larger households tended to have more contacts. Contact patterns differed significantly across time periods. Weekdays during the school semester showed the highest number of contacts, followed by weekdays during vacations, the Lunar New Year holidays, and weekends. During the Lunar New Year, there was an increase in contacts with extended family members, enhancing subnational social mixing. These findings provide valuable insights into social contact patterns in South Korea, which can improve our understanding of disease transmission and aid in developing region-specific epidemiological models.

The unpredictable timing of infectious disease outbreaks poses significant challenges for public health preparedness. For many pathogens, this unpredictability is due to uncertainty regarding introduction rates—the frequency with which the pathogen is introduced into at-risk populations. We present three model-driven advances toward quantifying pathogen introduction rates and their effects on outbreak timing and size. Our method relies on the assumption that pathogen introductions can only cause large outbreaks when population immunity is sufficiently low (i.e., the reproduction number R(t) > 1). First, we demonstrate that, for pathogens that cause lifelong immunity, introduction rates can be estimated from age-structured serological data. Second, we estimate annual rates of chikungunya virus (CHIKV, a mosquito-borne pathogen) introductions into 17 populations in Africa and Asia using serological data collected between 1973 and 2015. Our median estimates ranged from 1 to 70 CHIKV introductions per 10 million people per year. Finally, we used simulations to show how the introduction rate of a pathogen can shape its transmission patterns over time in affected populations. A lower introduction rate allows population immunity to wane between introductions, leading to large but infrequent outbreaks. In contrast, a higher introduction rate causes frequent low-level transmission, resulting in elevated population immunity that precludes large outbreaks. Together, these results illustrate how age-structured serology, a common type of epidemiological data, can be leveraged to better understand both historical and future transmission patterns in different populations.

As mpox spread continues across several endemic and non-endemic countries around the world, vaccination has become an integral part of the global response to control the epidemic. Some vaccines have been recommended for use against mpox by the World Health organization (WHO). As the roll-out of mpox vaccines continue across the globe, it is imperative to develop mathematical models to support public health officials and governments agencies in optimizing vaccination strategies to curtail the resurgence of mpox. In this article, we develop a compartmental mathematical model to investigate the impact of vaccination in controlling a potential mpox resurgence in Canada. The model categorizes individuals into high- and low-risk groups and incorporates pre-exposure vaccination in the high-risk group and post-exposure vaccination in the high- and low-risk groups. The vaccine-free version of the model was calibrated to the daily reported cases of mpox in Canada from April to October 2022, from which we estimated key model parameters, including the sexual and non-sexual transmission rates. Furthermore, we calibrated the full model to the daily reported cases of mpox in Canada in 2024, to estimate the current mpox vaccination rates in Canada. Our results highlight the importance of pre-exposure vaccination in the high-risk group on controlling a potential resurgence of mpox in Canada, and the minimal effects of post-exposure vaccination in the high- and low-risk groups on the outbreak. In addition, our model predicts the possibility of mpox becoming endemic in Canada, in the absence of pre-exposure vaccination in the high-risk group. Overall, our modeling result suggests that pre-exposure vaccination in the high-risk group is crucial in controlling mpox outbreak in Canada. Stepping up this vaccination is sufficient to avert a potential mpox resurgence in Canada.

A deterministic nonlinear mathematical model describing the population dynamics for Dengue and Chikungunya virus taken into consideration the effect of misdiagnosis due to the co-endemicity of the two viruses in the human population. It is necessary to understand the most important parameters involved in their dynamics that may help in developing strategies for prevention, control and joint treatments. The model is rigorously analyzed qualitatively and thresholds for eradication are established.

Mosquito-borne diseases, such as malaria, remain a major global health challenge, necessitating the exploration of innovative vector control strategies. Naturally occurring entomopathogenic fungi have been shown to reduce mosquito lifespan, but their slow-acting nature has limited their practical application. Advances in biotechnology have led to the development of transgenic fungus strains (this study will focus on Metarhizium pingshaense strain) engineered to express insecticidal toxins, significantly increasing their efficacy against malaria vector mosquitoes. To our knowledge, this is the first deterministic model designed to assess the impact of fungal-based mosquito control. The proposed model accounts for multiple transmission pathways of the fungal infection, including mating-based transmission from infected males to females and indirect transmission via contact with infectious mosquito carcasses. The model is analyzed to determine equilibrium states, local stability conditions, and the reproduction number. Numerical simulations explore various release scenarios, evaluating the effectiveness of periodic versus continuous fungal release in different ecological settings. The results indicate that transgenic Metarhizium pingshaense has the potential to significantly reduce mosquito populations, particularly when release strategies are optimized.

Human Papillomavirus (HPV) is a major sexually-transmitted infection that causes various cancers and genital warts in humans. In addition to accounting for about 99% of cervical cancer cases, it significantly contributes to anal, penile, vaginal, and head and neck cancers. Although HPV is vaccine-preventable (and highly efficacious vaccines exist for preventing infection with some of the most oncogenic HPV subtypes in the targeted population), the disease continues to cause major public health burden globally (largely due to inequity in access to the main control resources (i.e., access to Pap smear and vaccination) and low vaccination coverage in most communities that implement routine HPV vaccination). This lecture is based on the use of a new mathematical model (for the natural history of HPV, together with the associated neoplasia) for assessing the combined population-level impacts of Pap cytology screening and vaccination against the spread of HPV in a heterogeneous (heterosexual and homosexual) population. The model, which takes the form of a deterministic system of nonlinear differential equations, will be calibrated and validated using HPV-related cancer data from South Korea. Theoretical and numerical simulation results will be presented. Conditions for achieving vaccine-derived herd-immunity threshold (for achieving HPV elimination in Korea) will be derived.

A genome consists of a long string of four letters (bases A, T, C, G). How the information of biochemical processes stored in this string of bases is an open question. Are their higher order structures, such as, words, sentences, semantics, and a grammar in the DNA language (compositional complexity)? DNA from different species exhibit differences in global sequence composition, and this is used as markers to align larger sequences - grouping of genomes based on homology. Classification of genomes through similarity and dissimilarity is at the heart of Phylogenetics/Genomic Epidemiology. It uses several statistical-mathematical methods to align and compare the base sequences of multiple genomes, which are both computational resource intensive and time consuming for similar sequences. We develop and use an “alignment-free” method based on the Chaos-Game-Representation (CGR) of Statistical Physics, to successfully classify very closely related genomes of sub and sub-sub-species of HIV1 and mutants of Covid19. This useful approach requires less computational resources and time for analysis.

Biological systems evolve across different spatial and temporal scales. Modeling such complex systems gives rise to multi-scale differential equations that may be written as ODEs, PDEs, DDEs, SDEs, or Difference Equations. Particularly, vector-borne disease models are often described using ordinary differential equations with multiple time scales, which can involve singular perturbations—situations where rapid transitions or significant changes in system behavior occur due to small parameter variations or the interaction between fast and slow dynamics. To analyze these multi- time scale problems, we employ tools such as Geometric Singular Perturbation Theory (GSPT), Tikhonov’s Theorem, and Fenichel’s Theory. These methods provide insights into complex phenomena, including the loss of normal hyperbolicity and other intricate behaviors. Particularly, applying singular perturbation theory to vector-borne diseases allows us to reduce the dynamics to a one-time scale and understand their dynamics. To illustrate this, we present a Zika virus model and apply Tikhonov’s theorem and GSPT to investigate the model’s asymptotic behavior. Additionally, we conduct a bifurcation analysis to explore how the system’s behavior changes with variations in the parameter . We illustrate the various qualitative scenarios of the reduced system under singular perturbation. We will show that the fast–slow models, even though in nonstandard form, can be studied by means of GSPT.

Despite being curable and preventable, tuberculosis (TB) still causes the highest mortality rates in the human population. However, the number of TB cases significantly reduced globally in 2020, according to the Global Tuberculosis Report by the World Health Organization, coinciding with the COVID-19 pandemic. These reductions in TB cases are likely due to a complex interplay between disruptions in TB health services and the case counts resulting from COVID-19. We developed a compartmental model for the co-infection of tuberculosis and COVID-19 in the human population to assess the impact of medical resources, mobility, under-reporting, and the social behavior (follow social distancing and face mask) of infected individuals with either disease. We computed the basic reproduction numbers for TB alone, COVID-19 alone, and the co-infection scenario. Additionally, key parameters and basic reproduction numbers were estimated by utilizing case studies from low-income, middle-income, and high-income countries in a multi-patch scenario. Our results indicate that increased social behavior among infected individuals significantly reduces the number of co-infected individuals. The impact of mobility was assessed using a two-patch model with emigration and immigration rates. It shows that the mobility of unreported infectious individuals significantly increases both active cases of TB and COVID-19. This study provides significant recommendations for medical providers and public health officials regarding TB elimination in high-income countries and TB reduction in lower-income countries with a high disease burden. The findings are also relevant for studying TB in the context of future pandemic scenarios.

Chlamydia trachomatis (CT) infection is the most reported bacterial sexually transmitted infection (STI) in the United States (US). Despite many cases being asymptomatic, infection can lead to complications such as pelvic inflammatory disease (PID) in females, and infertility in both females and males. We developed an agent-based transmission model to evaluate the impact of a potential CT vaccine on the prevalence of CT infections and associated PID in the US population. The model simulates an evolving sexual network of 10,000 sexually active agents aged 15–54, including heterosexuals, female sex workers, and men who have sex with men, following Susceptible–Exposed–Infected–Recovered–Susceptible (SEIRS) transmission dynamics. A key strength of the model is its rigorous two-step calibration procedure, which first matches real CT prevalence by age and sex, followed by real PID prevalence by age in the US. This ensures realistic alignment with epidemiological patterns. The model incorporates both vaccination and test-and-treat strategies, enabling direct comparisons between interventions. We then evaluated the impact of different scenarios of vaccination coverage and targeting, assuming a vaccine with 80% efficacy against infection and a 5-year duration of protection. The results demonstrate a gender-neutral vaccine recommendation is projected to achieve the highest impact in reducing CT prevalence and PID burden, even with a moderate vaccination coverage. Beyond CT, this is flexible, computationally efficient framework is adaptable to study other STIs and assess the effectiveness of various intervention strategies, given appropriate epidemiological and behavioral data. By providing actionable insights, this framework serves as a decision-support tool for policymakers, public health officials, and vaccine developers.

A deterministic nonlinear mathematical model describing the population dynamics for Dengue and Chikungunya virus taken into consideration the effect of misdiagnosis due to the co-endemicity of the two viruses in the human population. It is necessary to understand the most important parameters involved in their dynamics that may help in developing strategies for prevention, control and joint treatments. The model is rigorously analyzed qualitatively and thresholds for eradication are established.

Mosquito-borne diseases, such as malaria, remain a major global health challenge, necessitating the exploration of innovative vector control strategies. Naturally occurring entomopathogenic fungi have been shown to reduce mosquito lifespan, but their slow-acting nature has limited their practical application. Advances in biotechnology have led to the development of transgenic fungus strains (this study will focus on Metarhizium pingshaense strain) engineered to express insecticidal toxins, significantly increasing their efficacy against malaria vector mosquitoes. To our knowledge, this is the first deterministic model designed to assess the impact of fungal-based mosquito control. The proposed model accounts for multiple transmission pathways of the fungal infection, including mating-based transmission from infected males to females and indirect transmission via contact with infectious mosquito carcasses. The model is analyzed to determine equilibrium states, local stability conditions, and the reproduction number. Numerical simulations explore various release scenarios, evaluating the effectiveness of periodic versus continuous fungal release in different ecological settings. The results indicate that transgenic Metarhizium pingshaense has the potential to significantly reduce mosquito populations, particularly when release strategies are optimized.

Human Papillomavirus (HPV) is a major sexually-transmitted infection that causes various cancers and genital warts in humans. In addition to accounting for about 99% of cervical cancer cases, it significantly contributes to anal, penile, vaginal, and head and neck cancers. Although HPV is vaccine-preventable (and highly efficacious vaccines exist for preventing infection with some of the most oncogenic HPV subtypes in the targeted population), the disease continues to cause major public health burden globally (largely due to inequity in access to the main control resources (i.e., access to Pap smear and vaccination) and low vaccination coverage in most communities that implement routine HPV vaccination). This lecture is based on the use of a new mathematical model (for the natural history of HPV, together with the associated neoplasia) for assessing the combined population-level impacts of Pap cytology screening and vaccination against the spread of HPV in a heterogeneous (heterosexual and homosexual) population. The model, which takes the form of a deterministic system of nonlinear differential equations, will be calibrated and validated using HPV-related cancer data from South Korea. Theoretical and numerical simulation results will be presented. Conditions for achieving vaccine-derived herd-immunity threshold (for achieving HPV elimination in Korea) will be derived.

A genome consists of a long string of four letters (bases A, T, C, G). How the information of biochemical processes stored in this string of bases is an open question. Are their higher order structures, such as, words, sentences, semantics, and a grammar in the DNA language (compositional complexity)? DNA from different species exhibit differences in global sequence composition, and this is used as markers to align larger sequences - grouping of genomes based on homology. Classification of genomes through similarity and dissimilarity is at the heart of Phylogenetics/Genomic Epidemiology. It uses several statistical-mathematical methods to align and compare the base sequences of multiple genomes, which are both computational resource intensive and time consuming for similar sequences. We develop and use an “alignment-free” method based on the Chaos-Game-Representation (CGR) of Statistical Physics, to successfully classify very closely related genomes of sub and sub-sub-species of HIV1 and mutants of Covid19. This useful approach requires less computational resources and time for analysis.

Biological systems evolve across different spatial and temporal scales. Modeling such complex systems gives rise to multi-scale differential equations that may be written as ODEs, PDEs, DDEs, SDEs, or Difference Equations. Particularly, vector-borne disease models are often described using ordinary differential equations with multiple time scales, which can involve singular perturbations—situations where rapid transitions or significant changes in system behavior occur due to small parameter variations or the interaction between fast and slow dynamics. To analyze these multi- time scale problems, we employ tools such as Geometric Singular Perturbation Theory (GSPT), Tikhonov’s Theorem, and Fenichel’s Theory. These methods provide insights into complex phenomena, including the loss of normal hyperbolicity and other intricate behaviors. Particularly, applying singular perturbation theory to vector-borne diseases allows us to reduce the dynamics to a one-time scale and understand their dynamics. To illustrate this, we present a Zika virus model and apply Tikhonov’s theorem and GSPT to investigate the model’s asymptotic behavior. Additionally, we conduct a bifurcation analysis to explore how the system’s behavior changes with variations in the parameter . We illustrate the various qualitative scenarios of the reduced system under singular perturbation. We will show that the fast–slow models, even though in nonstandard form, can be studied by means of GSPT.

Despite being curable and preventable, tuberculosis (TB) still causes the highest mortality rates in the human population. However, the number of TB cases significantly reduced globally in 2020, according to the Global Tuberculosis Report by the World Health Organization, coinciding with the COVID-19 pandemic. These reductions in TB cases are likely due to a complex interplay between disruptions in TB health services and the case counts resulting from COVID-19. We developed a compartmental model for the co-infection of tuberculosis and COVID-19 in the human population to assess the impact of medical resources, mobility, under-reporting, and the social behavior (follow social distancing and face mask) of infected individuals with either disease. We computed the basic reproduction numbers for TB alone, COVID-19 alone, and the co-infection scenario. Additionally, key parameters and basic reproduction numbers were estimated by utilizing case studies from low-income, middle-income, and high-income countries in a multi-patch scenario. Our results indicate that increased social behavior among infected individuals significantly reduces the number of co-infected individuals. The impact of mobility was assessed using a two-patch model with emigration and immigration rates. It shows that the mobility of unreported infectious individuals significantly increases both active cases of TB and COVID-19. This study provides significant recommendations for medical providers and public health officials regarding TB elimination in high-income countries and TB reduction in lower-income countries with a high disease burden. The findings are also relevant for studying TB in the context of future pandemic scenarios.

Chlamydia trachomatis (CT) infection is the most reported bacterial sexually transmitted infection (STI) in the United States (US). Despite many cases being asymptomatic, infection can lead to complications such as pelvic inflammatory disease (PID) in females, and infertility in both females and males. We developed an agent-based transmission model to evaluate the impact of a potential CT vaccine on the prevalence of CT infections and associated PID in the US population. The model simulates an evolving sexual network of 10,000 sexually active agents aged 15–54, including heterosexuals, female sex workers, and men who have sex with men, following Susceptible–Exposed–Infected–Recovered–Susceptible (SEIRS) transmission dynamics. A key strength of the model is its rigorous two-step calibration procedure, which first matches real CT prevalence by age and sex, followed by real PID prevalence by age in the US. This ensures realistic alignment with epidemiological patterns. The model incorporates both vaccination and test-and-treat strategies, enabling direct comparisons between interventions. We then evaluated the impact of different scenarios of vaccination coverage and targeting, assuming a vaccine with 80% efficacy against infection and a 5-year duration of protection. The results demonstrate a gender-neutral vaccine recommendation is projected to achieve the highest impact in reducing CT prevalence and PID burden, even with a moderate vaccination coverage. Beyond CT, this is flexible, computationally efficient framework is adaptable to study other STIs and assess the effectiveness of various intervention strategies, given appropriate epidemiological and behavioral data. By providing actionable insights, this framework serves as a decision-support tool for policymakers, public health officials, and vaccine developers.

Deep learning has emerged in many fields in recent times where neural networks are used to learn and understand data. This study combines deep learning frameworks with epidemiological models and is aimed specifically at the creation and testing of DINNs with a view to modeling the infection dynamics of epidemics. Our research thus trains the DINN on synthetic data derived from an SI-SIR model designed for Avian influenza and shows the model’s accuracy in predicting extinction and persistence conditions. In the method, a twelve hidden layer model was constructed with sixty-four neurons per layer and ReLU activation function was used. The network is trained to predict the time evolution of five state variables for birds and humans over 50,000 epochs. The overall loss minimized to 0.000006, characterized by a combination of data and physics losses, enabling the DINN to follow the differential equations describing the disease progression.

During the early stages of the COVID-19 pandemic, it was important to minimize both medical and economic costs. In this study, we introduce a machine learning-based multi-objective optimization framework that can propose cost-effective social distancing strategies. Our approach finds Pareto solutions that balance different goals, like reducing infections and minimizing social distancing costs. Then, the cost-benefit analysis can adjust each cost factor—for example, value of statistical life (VSL), fatality rate, or GDP. We also provide an interactive web dashboard so that policymakers and the public can test various scenarios easily. We tested this framework on the COVID-19 pandemic in Korea. The results show that the difference between the cost-optimal strategy and implemented strategy is 10% in cost. Notably, our results reveal two distinct patterns in cost-optimal solutions. When social distancing cost is proportional to intervention intensity, an on-off lockdown strategy proves most economical. In contrast, when the cost increases sharply with intensifying social distancing, maintaining a consistently moderate level of intervention minimizes overall expenses. By letting people explore different cost settings and intervention strategies, this tool can support more balanced decisions during emerging infectious disease crises in the future.

The COVID-19 pandemic took its toll not only on the physical health of those who lived through it, but also on their mental health. I will present a model of the direct and indirect effects of COVID-19 and the associated public policies on mental health. This is an SIRS model of COVID-19, with compartments for mild, acute, and chronic COVID-19 infections and additional compartments for populations with mental health symptoms. Parameters are determined based on fitting from the first wave of COVID-19 in the New York state population, which includes several changes in local government policy, e.g, lockdown orders, which have an effect on the rate at which mental health systems develop. Finally, an additional “delay” term is included in the model to account for the delay between lockdowns going into effect and individuals developing mental health symptoms. The goal of our analysis is to understand how government policy in response to a pandemic can seek to maximize the population's quality-adjusted life years (QALY), which is a measure not only of lifespan, but also the quality of the years lived. I will present some preliminary results suggesting the optimal timing and strength of government lockdown mandates.

After an individual has been infected by a pathogen, T lymphocytes store information about the pathogen. Consequently, upon reinfection by the same pathogen, an immune response memory is triggered. This immune memory allows the body to react very quickly against the pathogen. Indeed, when an individual recovers from a virus, sometimes the individual acquires full immunity. In some cases, the individual's immunity persists for some period, after which it decreases progressively and can even disappear. I will then present a stochastic epidemic model with memory on the previous infection, incorporating varying infectivity and waning immunity. More precisely, we will present a functional law of large numbers when the size of the population tend to infinity. We will also present results on the behaviour of the epidemic, more precisely the threshold for the existence of an endemic equilibrium, and study the stability of the endemic equilibrium.

Symptom propagation occurs when an individual’s symptom severity is correlated with the symptom severity of the individual who infected them. Determining whether - and to what extent - these correlations exist requires data-driven methods. In this study, we use synthetic data to determine the types of data required to estimate the strength of symptom propagation and investigate the effect of reporting bias on these estimates. We found that even a relatively small number of contact tracing data points was sufficient to gain a reasonable estimate for the strength of symptom propagation. Increasing the number of contact tracing data points further improved our estimates. In contrast, population incidence alone was insufficient to accurately estimate the symptom propagation parameters, even with a large number of data points. Nonetheless, concurrently using population incidence data with contact tracing data led to increased accuracy when estimating the overall disease severity. We then considered the effect of severe cases being more likely to be reported in the contact tracing data. When contact tracing data alone was used, we found that our estimates for the strength of symptom propagation were robust to all reporting bias scenarios considered. However, the reporting bias led us to overestimate the overall disease severity. Using population incidence data in addition to contact tracing data reduced the error in disease severity but at the cost of increasing the error in the strength of symptom propagation when reporting bias was in both primary and secondary cases. Consequently, these errors led to us sometimes finding support for symptom propagation, even when the synthetic data was generated without.

We aim to compare results for the predicted reduction in vectorial capacity caused by pyrethroid and pyrethroid-piperonyl butoxide insecticide treated nets (ITNs) between semi-field Ifakara Ambiant Chamber tests (I-ACT) and experimental hut experiments. Mathematical modelling and Bayesian inference frameworks estimated ITN effects on mosquito behavioural endpoints (repelled, killed before/after feeding) to predict reductions in Anopheles gambiae’s vectorial capacity for Plasmodium falciparum transmission. The reduction in biting estimates are generally greater for I-ACT, possibly due to lower mosquito aggression: Although I-ACT vectors are probing before release, experimental hut vectors are actively seeking a blood meal. I-ACT estimates higher probability of killing vectors which have fed, while experimental huts show greater killing before feeding, possibly due to their open-system design, where vectors can contact the net, then attempt to exit and get trapped. This is supported by most of the mosquitoes being caught before feeding being in the exit trap. While the I-ACT is a closed system, were vectors cannot exit or be trapped, increasing the likelihood of returning to host-seeking and feeding. Despite these differences, both methods yielded similar predictions for the overall reduction in vectorial capacity. Results suggest that I-ACT provides a good initial assessment of the impact of adulticide modes of action of these nets. Challenges of semi-field experiments include how to model the change in efficacy from practical use over time. However, important advantages include the ability to easily trial different strains of vector (including different resistance levels) and allowing rapid data collection. Parameterising models with location-specific bionomic parameters allows for setting -specific predictions of the impact of different nets, with the potential to include additional modes of action for other active ingredients.

We aim to compare results for the predicted reduction in vectorial capacity caused by pyrethroid and pyrethroid-piperonyl butoxide insecticide treated nets (ITNs) between semi-field Ifakara Ambiant Chamber tests (I-ACT) and experimental hut experiments. Mathematical modelling and Bayesian inference frameworks estimated ITN effects on mosquito behavioural endpoints (repelled, killed before/after feeding) to predict reductions in Anopheles gambiae’s vectorial capacity for Plasmodium falciparum transmission. The reduction in biting estimates are generally greater for I-ACT, possibly due to lower mosquito aggression: Although I-ACT vectors are probing before release, experimental hut vectors are actively seeking a blood meal. I-ACT estimates higher probability of killing vectors which have fed, while experimental huts show greater killing before feeding, possibly due to their open-system design, where vectors can contact the net, then attempt to exit and get trapped. This is supported by most of the mosquitoes being caught before feeding being in the exit trap. While the I-ACT is a closed system, were vectors cannot exit or be trapped, increasing the likelihood of returning to host-seeking and feeding. Despite these differences, both methods yielded similar predictions for the overall reduction in vectorial capacity. Results suggest that I-ACT provides a good initial assessment of the impact of adulticide modes of action of these nets. Challenges of semi-field experiments include how to model the change in efficacy from practical use over time. However, important advantages include the ability to easily trial different strains of vector (including different resistance levels) and allowing rapid data collection. Parameterising models with location-specific bionomic parameters allows for setting -specific predictions of the impact of different nets, with the potential to include additional modes of action for other active ingredients.

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

Deterministic and stochastic models, though often used to describe the same biological systems, can yield qualitatively different predictions. In particular, deterministic bistability does not necessarily imply stochastic bimodality, and vice versa. Multistability and multimodality are typically seen as indicators of distinct system behaviors, often inferred from stochastic simulation trajectories. In this talk, I explore the disconnect between probability modes and behavioral modes in the context of transient bimodality—also known as “adiabatic explosions”—which refers to metastable probability modes that do not correspond to distinct dynamical behaviors at the trajectory level. I present recent findings that link the emergence of these transient modes to a breakdown of the central limit theorem, specifically in the context of first-passage time distributions to absorbing states. I conclude by discussing how transient bimodality challenges conventional interpretations of system behavior in biological contexts and highlight conditions under which this phenomenon becomes particularly relevant.

Mathematical models of ecological and epidemiological systems often focus on asymptotic dynamics, such as equilibria and periodic orbits. However, many systems exhibit long transient behaviors where certain variables of interest remain in a slowly evolving state for an extended period before undergoing rapid change. These transient dynamics can have significant implications for population persistence, disease outbreaks, and ecosystem stability. In this work, we analyze long-lasting and slowly varying transient dynamics in discrete-time systems. We extend previous theoretical frameworks by identifying conditions under which an observable of the system can exhibit prolonged transients and derive criteria for characterizing these dynamics. Our results show that specific points in the state space, analogous to transient centers in continuous-time systems, can generate and sustain long transients in discrete-time models. We further demonstrate how these properties manifest in predator-prey models and epidemiological systems, particularly in contexts where populations or disease prevalence remain low for an extended period before experiencing a sudden shift. These findings provide a foundation for understanding and predicting long transients in discrete-time ecological and epidemiological models. (Joint Work with FMG Magpantay)

Reinforcement learning (RL) is a powerful machine learning paradigm in which an RL agent learns to discover optimal strategies in uncertain environments. The RL control strategy has achieved remarkable success in complex tasks such as playing Chess, Go, and StarCraft. For RL, the prevailing application interface (API) standard is Gymnasium, a Python library [1]. Agent-based (AB) modeling is a mathematical, dynamical system modeling approach where the parts of the system, the so-called agents, autonomously act according to agent-type specific rules. PhysiCell is an AB modeling framework written in C++ and was implemented to model multicellular systems based on Newtonian physics. Cells are the agents. The cell type specifies the rule set the agents apply. Tissue structure emerges from the cell interactions. Substrates like oxygen can be modeled with the integrated BioFVM diffusive transport solver. Additionally, intracellular models can be integrated into cell agents [2]. The resulting AB models are 2 or 3-dimensional, off-lattice, center-based, and multiscale in space and time. In this talk, we will introduce PhysiGym, a well-documented and on all major operating systems tested open-source framework written in C++ and Python that allows to control PhysiCell models over the Gymnasium API. After a brief introduction to AB models and RL, we will discuss the implementation and obtained results from our tumor microenvironment model and the RL algorithms we applied to the model. In the future, PhysiGym can be used to learn from simulations possible mechanisms that might explain how biology systems react to similar real-world control. Furthermore, if cancer patient digital twins are written as PhysiCell models, PhysiGym could ultimately be used by oncologists to explore RL reward functions to improve treatment efficacy, reduce side effects, and slow or prevent resistance. References: [1] https://gymnasium.farama.org/ , [2] https://PhysiCell.org

Chemical Reaction Networks (CRNs) provide a powerful framework for modeling interactions between multiple chemical species within complex biological pathways. Their dynamics can be described by the mass-action law, representing variations over time in species concentrations through large systems of ODEs. The study of CRNs modeling signaling mechanisms within single cells has been successfully applied to provide insight into oncogenic pathways, enabling more predictive cancer models and improved therapeutic strategies [1,2,3]. Developing a CRN requires selecting the involved species, defining their interactions, and identifying key parameters such as microscopic reaction rates. This talk presents a method to retrieve these rates, particularly in the context of enzyme kinetics. We consider a simple model, where an enzyme E binds to a substrate S, forming a complex C that dissociates into a product P while regenerating E. This process is governed by four ODEs with three reaction rates: forward k_f, reverse k_r, and catalytic k_cat, typically hard to determine experimentally. A common simplification leads to the Michaelis-Menten (MM) model, where enzyme kinetics is characterized by a single ODE depending on two measurable parameters: the Michaelis constant K_M, and the maximum reaction rate V_max. These parameters may be expressed as functions of microscopic rates and are more accessible experimentally. This talk addresses the inverse problem of estimating k_f and k_r from K_M and V_max. A computational algorithm for solving this problem is presented and analyzed, along with an estimate of the reconstruction accuracy, and numerical simulations demonstrating its potential for refining kinetic models in biological and biomedical research. [1] Sommariva et al., J. Math. Biol., 82 (6): 55, 2021. [2] Berra et al., J. Optim. Theory Appl., 200 (1): 404-427, 2024. [3] Sommariva et al., Front. Syst. Biol., 3: 1207898, 2023.

Most biological systems have underlying multiple spatial or temporal scales that require reduced-order models to capture their essential dynamics and analyze them. However, traditional model reduction techniques, such as Computational Singular Perturbation (CSP), rely on the availability of the governing or dynamical equations, which are often unknown from data in biomedical applications. To address this limitation, we propose a data-driven CSP framework that integrates Sparse Identification of Nonlinear Dynamics (SINDy) and Neural Networks to extract time-scale separated models directly from data. Our approach is validated on the Michaelis-Menten enzyme kinetics model, a well-established multiscale system, by identifying reduced models for the standard Quasi-Steady-State Approximation (sQSSA) and reverse Quasi-Steady-State Approximation (rQSSA). When the full model cannot be identified by SINDy due to noise, we use Neural Networks to estimate the Jacobian matrix, allowing CSP to determine the regions where reduced models are valid. We further analyze Partial Equilibrium Approximation (PEA) case, where the dynamics span both sQSSA and rQSSA regimes, requiring dataset splitting to accurately identify region-specific models. The results demonstrate that SINDy struggles in the presence of noise to identify full model from data that have underlying timescale evolution, but remains effective for identifying reduced models when dataset are partitioned correctly.

Although chemical reaction networks (CRN) provide performant and biophysically plausible models for explaining single-cell genomic data, inference of reaction network parameters in this setting usually assumes available data points can be viewed as independent samples from a steady state distribution. Less is known about how to perform efficient parameter inference in the case that there is a continuous-time data stream, which adds complexity like nontrivial correlations between samples from different times. In the continuous-time setting, one has two natural questions: (i) given a set of reactions that could plausibly explain the observed data stream, what are reasonable estimates of the associated reaction rate parameters? and (ii) what is the minimal set of reactions necessary to explain the data? Both questions can be formalized as Bayesian inference problems, with the former concerning the inference of a model-dependent parameter posterior, and the latter concerning ‘structure’ inference. If one can assume each possible reaction has a different stoichiometry vector, there is a well-known analytic solution to both problems; if reactions can have the same stoichiometry vector (i.e., there is reaction degeneracy), both problems become substantially more difficult, and no analytic solution is known. We present the first approximate analytic solution to both problems, which is valid when the number of observations becomes sufficiently large. In its regime of validity, this solution allows its user to avoid expensive likelihood computations that can involve summing over an exponentially large number of terms. We discuss interesting consequences of this solution, like the fact that ‘simpler’ models with fewer reactions are preferred over more complex ones, and the fact that the parameter posteriors of non-identifiable models are strongly prior-dependent.

Real-world systems, especially in biology, exhibit significant complexity and inherent limitations in observability. What methods can enhance our understanding of the mechanisms underlying their functionality? Additionally, how can we develop and test explanatory models within a stochastic environment? Evaluating the effectiveness of these models requires quantitative measurements of the disparity between model outputs and observed data. While mean-field, deterministic models have well-established approaches for such assessments, stochastic systems—particularly those constrained by multiple data types—need carefully designed quantitative comparison methods. Methods for inferring the parameters of stochastic models generally require analytical forms of the model solutions, large data sets, summary statistics, or assumptions on the distribution of model outputs. These approaches can be limiting if you wish to preserve the information in the variability of the data but you do not have sufficient data to reliably fit distributions or determine robust statistics. We present a hierarchical approach to develop a distance measure for the direct comparison of model output distributions to experimentally observed distributions, avoiding any assumptions about distributions and the need to choose summary statistics. Our distance measure allows for constraining the model with multiple experiments, not necessarily of the same type, such that each experiment constrains some, or all, of the model parameters. We use this distance for parameter estimation with our queuing model of intracellular GLUT4 translocation. We will explore some practical considerations when using the distance for parameter inference, such as the effects of model output sampling and experimental error. Fitting the queuing model to data allowed us to uncover a possible mechanism of GLUT4 sequestration and release in response to insulin. Authors: Brock D. Sherlock and Adelle C.F. Coster

Collective behaviour in animal groups are coordinated movements and interactions among members that aim to achieve a common goal. Whether these goals are for allocation of resources or defence from predators, the collective behaviour appears to be largely a group activity initiated by a member, known as the leader. In the absence of high-resolution spatio-temporal data, various qualitative studies offer a glimpse of how leader-follower interactions take place. For example, Nagy, et al., studied the average delay in response when pigeons change the direction inflight. Next, Bourjade, et al., studied the first mover and the succeeding order of movements of Przewalski's Horses. Furthermore, various studies on the collective motion in the animal kingdom offer mathematical models and infer how the interactions and decision making take place. Important questions arise during an event of coordinated motion in animals. During such an event, do individuals move according to a certain set of natural rules? Or certain patterns form due to the influence of a leader? How is this influence measured? Finally, how is influence transferred to other members of the group? In this study, we discuss the role of information theory to quantitatively uncover leader-follower relationship in a horse group. Specifically, we introduce concepts from information theory, specifically global and local transfer entropy being applied to a harem of horses. We will discuss their definitions, and how these key concepts are used to support causation in events. We will then discuss some important implications on how this technique can be used to analyse collective motion where data is scarce.

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

Alzheimer’s Disease (AD) is a progressive neurodegenerative disorder affecting approximately 10% of Americans over age 65, leading to memory loss, cognitive decline, and impaired daily function. Disease progression correlates with the spread of tau and amyloid-beta proteins, which aggregate into neurofibrillary tangles. While macroscopic whole-brain network models predict large-scale protein deposition patterns, they lack the specificity to capture individual disease progression. Conversely, microscale neuron-neuron models offer highly detailed biochemical aggregation and transport simulations but are computationally prohibitive for whole-brain parameter inference. In this work, we explore using machine learning to accelerate whole-brain simulations by approximating explicit solutions to the microscopic Two-Neuron Transport Model (TNTM), a partial differential equation describing tau flux along a neuron, incorporating biochemical aggregation, fragmentation, and transport. We simulate a single-edge model across physiological ranges of biochemical parameters and boundary conditions, then compare regression methods with varying levels of interpretability, from neural networks (low) to symbolic regression via PySR (high). Neural networks achieve the lowest error but lack biological insight. Linear and polynomial regression compute rapidly but yield high errors with limited interpretability. Symbolic regression achieves a balance between accuracy and transparency. This work demonstrates the potential of machine learning for computationally scalable AD modeling, opening avenues for patient-specific parameterization using AD data repositories.

We introduce a method to identify phase equations for heterogeneous oscillators beyond the weak coupling regime. This strategy is an extension of the theory from [Y. Park and D. Wilson, SIAM J. Appl. Dyn. Syst., 20 (2021), pp. 1464--1484] and yields coupling functions for N general limit-cycle oscillators with arbitrary types of coupling, with similar benefits as the classic theory of weakly coupled oscillators. These coupling functions enable the study of oscillator networks in terms of phase-locked states, whose stability can be determined using straightforward linear stability arguments. We demonstrate the utility of this approach by reducing and analyzing conductance-based thalamic neuron model. The reduction correctly predicts the emergence of new phase-locked states as a function of coupling strength and heterogeneity. We conclude with a brief remark on recent extensions to n:m phase-locking and N-body interactions.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Microinvasions in ductal carcinoma in situ (DCIS) are malignant cells that have broken through the basement membrane (BM) and extend into the stroma with no focus larger than 1 mm. Since microinvasions constitute the first step in the metastatic cascade, identifying the causes of microinvasions will help distinguish between progressors or non-progressors among the DCIS patients, thus improving treatment. The mechanical tumor-stroma interactions play an important role in this process. Studies have shown that elevated collagen stiffening, deposition, and fibril crosslinking are correlated with tumor aggressiveness and invasion in breast cancer. Therefore, here we present SilicoDCIS, a 2D off-lattice center-based agent-based model (ABM) of ductal carcinoma in situ (DCIS) growth and its interaction with the extracellular matrix (ECM) to investigate the mechanical conditions that may lead to tumor microinvasions. SilicoDCIS simulates the division, growth, and migration of tumor cells in DCIS while interacting with other cell types and the ECM. This includes the BM, the myoepithelial and epithelial cell layers, and the collagen in the ECM. The ECM was modeled as a vector field, where the direction of each vector gives the orientation of a collagen bundle, and the vector magnitude is related to the bundle density. The growing DCIS can remodel the ECM (density and orientation), and in turn, the ECM applies a reciprocal force (proportional to the local collagen density) opposite to the tumor growth. With SilicoDCIS, we studied the mechanical effects of cancer cell proliferation and migration on the BM and the ECM. We found that higher cell migration force leads to increased BM stress and ECM density (on the tumor edges where cells migrate) and that the escape of the migrating cells from the duct vs. their intraductal confinement depends on cell speed. SilicoDCIS may provide insights into the mechanics of DCIS microinvasions to guide the design of future experiments.

Compound birth-death processes are widely used to model the age-incidence curves of many cancers [1]. There are efficient schemes for directly computing the relevant probability distributions in the context of linear multi-stage clonal expansion (MSCE) models [2]. However, these schemes have not been generalised to models on arbitrary graphs, forcing the use of either full stochastic simulations or mean-field approximations, which can become inaccurate at late times or old ages [3, 4]. Here, we present a numerical integration scheme for directly computing survival probabilities of a first-order birth-death process on an arbitrary directed graph, without the use of stochastic simulations. As a concrete application, we show that this new numerical method can be used to infer the parameters of an example graphical model from simulated data.

Melanoma is the most aggressive type of skin cancer, yet survival rates are excellent if it is diagnosed early. However, if metastasis occurs, five-year survival rates drop significantly. During the early stages of tumour initiation, melanoma cells form clusters within the primary tumour which promote metastasis. In the absence of biological tools to visualise cluster formation at primary tumour sites, we develop mathematical models to generate mechanistic insight into their formation. For this work we utilise in vitro data for two distinct melanoma cell phenotypes, one more proliferative and the other more invasive. This data consists of experiments for each phenotype individually, resulting in homogeneous clusters, as well as mixtures of the two phenotypes, resulting in heterogeneous clusters. We develop a series of differential-equation-based models using a coagulation-fragmentation-proliferation framework to describe the growth dynamics of homogeneous clusters, incorporating different functional forms for cell proliferation and cluster splitting. We then extend these models to describe the formation of heterogeneous cell clusters by considering both cluster size and phenotypic composition. We fit the models to experimental data, using a Bayesian framework to perform parameter inference and information criteria to perform model selection. In this way, we characterise and quantify differences in the clustering behaviour of two melanoma phenotypes in homogeneous and heterogeneous clusters, particularly the cluster coagulation, proliferation, and splitting rates. We find that the coagulation rate for the invasive phenotype is much larger than that for the proliferative phenotype, and evaluate how well different modelling assumptions fit the data in order to increase our understanding of the mechanisms driving metastasis. In future work, the models will be used to inform further experiments and, in particular, to suggest and test strategies for inhibiting metastasis.

A hallmark of colorectal cancer is the structural disruption of the colonic tissue, a process correlated with disease progression. Intestinal crypts, glands essential for homeostasis, lose their tubular morphology - and function - due to uncontrolled cell proliferation and tissue invasion. Evaluating this deterioration in biopsied samples is critical for both patient diagnosis and prognosis. Histopathological methods are essential for assessing colorectal cancer status, but their precision and reproducibility can be improved. Spatial biology provides a mathematical framework to analyse the structural properties of biological data; in this work, we apply techniques from topological data analysis and network science to quantify architectural changes in colorectal cancer progression. Using cell point clouds derived from immunohistochemistry imaging, we construct cell networks that encode topological tissue features. We employ these networks to segment large, imaged samples into smaller, biologically meaningful regions of interest that preserve tissue architecture. We compare the performance of our approach to conventional segmentation methods such as quadrat division. Within these segmented regions, we further employ methods from persistent homology to quantify tissue structure, with the long-term goal of identifying novel biomarkers of disease progression.

Cell-cell fusion is a known phenomenon throughout the human body. It characterizes a wide range of physiological and pathological processes, ranging from placentation and embryogenesis to cancer stem cell (CSC) formation. There is increasing evidence that cell-cell fusion can play key roles in the development and progression of cancer, particularly by increasing intratumor heterogeneity and potentiating somatic evolution. There are many unanswered questions surrounding the characteristics that define cancer cell-cell fusion events, their frequency in in vivo tumor conditions, and whether or not cell-cell fusion is a universal phenomenon across cancer. Using a combination of in vitro and in silico approaches, we can begin to answer some of these questions. We have developed a preliminary cellular automata model using HAL to evaluate the effect of variable cell-cell fusion rates and behaviors under a range of tumor microenvironmental conditions. By comparing our spatial model to a suite of ordinary differential equations, we can begin to estimate the effects of cell-cell fusion on the genomic heterogeneity and malignancy potential of cancers in vivo. I demonstrate the importance of improving fusion rate estimates using the simplest iteration of an in silico cellular automata model (coined SimpleFusion). The preliminary SimpleFusion model results illustrate how much the impact of cell fusion, as measured by the percentage of cells that have had a fusion event in their lineage, changes between orders of magnitude of fusion rates. Corresponding ODE models demonstrate similar results despite the lack of encoded spatial information. By studying these two types of models (ABM, ODEs) in combination, we can begin to understand what parameters most directly define the cell-cell fusion population dynamics in our in vitro fusion experiments and, in turn, in vivo conditions as well.

Blood cancers pose a growing medical and economic challenge in aging societies. Every day, the human bone marrow (BM) generates more than 100 billion blood cells. This process is driven by hematopoietic stem cells (HSCs), which retain their ability to proliferate and self-renew throughout life. However, over time, HSCs accumulate mutations that may lead to malignant transformation, as seen in acute myeloid leukemia (AML), one of the most aggressive cancers. Even in healthy individuals, the BM undergoes age-related changes, including a decline in cell numbers, remodeling of the BM micro-environment, and a bias in HSC differentiation. Emerging evidence suggests that these alterations create a favorable environment for the expansion of mutated cells, thereby promoting blood cancer development and progression. Mathematical and computational models facilitate our understanding of how BM aging contributes to malignant cell growth. We propose nonlinear ordinary differential equation models to describe blood cell formation and clonal competition in the human BM. The models incorporate micro-environmental and systemic feedback loops and are informed by data from both healthy individuals and cancer patients. Our findings suggest that the age-related decline in HSC self-renewal, combined with increased chronic inflammation (inflammaging), makes the BM more susceptible to the expansion of mutated cells and at the same time impairs treatment response. Through mathematical analysis, quantitative simulations, and patient data fitting, we study the following questions: 1. How do HSC proliferation & self-renewal change during physiological aging? 2. How do age-related alterations in healthy BM contribute to blood cancer development? 3. What is the impact of chronic inflammation on HSC function and blood cancer progression? 4. How do age-related BM changes affect treatment responses, e.g., in AML patients? 5. How could treatment protocols be adapted to elderly patients?

Glioblastoma multiforme (GBM) is a fast-growing and deadly brain tumor due to its ability to aggressively invade the nearby brain tissue. A host of mathematical models in the form of reaction–diffusion equations have been formulated and studied in order to assist clinical assessment of GBM growth and its treatment prediction. To better understand the speed of GBM growth and form, we propose a two population reaction–diffusion GBM model based on the ‘go or grow’ hypothesis. Our model is validated by in vitro data and assumes that tumor cells are more likely to leave and search for better locations when resources are more limited at their current positions. Our findings indicate that the tumor progresses slower than the simpler Fisher model, which is known to overestimate GBM progression. Moreover, we obtain accurate estimations of the traveling wave solution profiles under several plausible GBM cell switching scenarios by applying the approximation method introduced by Canosa.

Patients with inflammatory bowel disease (IBD) face an elevated risk of colorectal cancer (CRC), necessitating lifelong surveillance to find and remove precancers before they become malignant. Current one-size-fits-all approaches are inadequate and tailored strategies that consider cancer evolution are needed. To address this, we developed a mechanistic framework of IBD-CRC progression. Our multi-type branching process model accounts for IBD onset, mutational processes, and both precancerous (adenoma/dysplasia) and malignant clonal expansion. Initial parameter estimation for mutation and growth rates when fitting the multi-stage clonal expansion model to epidemiological IBD-CRC data yielded similar estimates to those found previously in sporadic CRC but suggest higher mutation rates and slightly lower growth rates in IBD. However, this data may not perfectly represent the natural history, as surveillance colonoscopy with lesion removal and colectomy alter the observable progression. Further, fitting to cancer incidence data alone presents parameter identifiability issues, restricting our initial fit to four parameters. To address these limitations, our study draws upon extensive clinical data from the U.S. Veterans Health Administration, employing validated methods using large language models to construct high-quality datasets with detailed information on surveillance colonoscopy timing, colectomies, and intermediate lesions extracted from pathology reports. To integrate these data, we developed a complementary fast simulation model, which will be released as an R package. This simulation model incorporates clinical interventions, such as colonoscopy with size-dependent lesion removal. Our combined analytical and simulation approach captures the complex precancerous evolution in IBD, providing a quantitative foundation for more effective, personalized surveillance guidelines. Further, this approach can be adapted to improve surveillance in the general population.

Blood cancers pose a growing medical and economic challenge in aging societies. Every day, the human bone marrow (BM) generates more than 100 billion blood cells. This process is driven by hematopoietic stem cells (HSCs), which retain their ability to proliferate and self-renew throughout life. However, over time, HSCs accumulate mutations that may lead to malignant transformation, as seen in acute myeloid leukemia (AML), one of the most aggressive cancers. Even in healthy individuals, the BM undergoes age-related changes, including a decline in cell numbers, remodeling of the BM micro-environment, and a bias in HSC differentiation. Emerging evidence suggests that these alterations create a favorable environment for the expansion of mutated cells, thereby promoting blood cancer development and progression. Mathematical and computational models facilitate our understanding of how BM aging contributes to malignant cell growth. We propose nonlinear ordinary differential equation models to describe blood cell formation and clonal competition in the human BM. The models incorporate micro-environmental and systemic feedback loops and are informed by data from both healthy individuals and cancer patients. Our findings suggest that the age-related decline in HSC self-renewal, combined with increased chronic inflammation (inflammaging), makes the BM more susceptible to the expansion of mutated cells and at the same time impairs treatment response. Through mathematical analysis, quantitative simulations, and patient data fitting, we study the following questions: 1. How do HSC proliferation & self-renewal change during physiological aging? 2. How do age-related alterations in healthy BM contribute to blood cancer development? 3. What is the impact of chronic inflammation on HSC function and blood cancer progression? 4. How do age-related BM changes affect treatment responses, e.g., in AML patients? 5. How could treatment protocols be adapted to elderly patients?

Glioblastoma multiforme (GBM) is a fast-growing and deadly brain tumor due to its ability to aggressively invade the nearby brain tissue. A host of mathematical models in the form of reaction–diffusion equations have been formulated and studied in order to assist clinical assessment of GBM growth and its treatment prediction. To better understand the speed of GBM growth and form, we propose a two population reaction–diffusion GBM model based on the ‘go or grow’ hypothesis. Our model is validated by in vitro data and assumes that tumor cells are more likely to leave and search for better locations when resources are more limited at their current positions. Our findings indicate that the tumor progresses slower than the simpler Fisher model, which is known to overestimate GBM progression. Moreover, we obtain accurate estimations of the traveling wave solution profiles under several plausible GBM cell switching scenarios by applying the approximation method introduced by Canosa.

Patients with inflammatory bowel disease (IBD) face an elevated risk of colorectal cancer (CRC), necessitating lifelong surveillance to find and remove precancers before they become malignant. Current one-size-fits-all approaches are inadequate and tailored strategies that consider cancer evolution are needed. To address this, we developed a mechanistic framework of IBD-CRC progression. Our multi-type branching process model accounts for IBD onset, mutational processes, and both precancerous (adenoma/dysplasia) and malignant clonal expansion. Initial parameter estimation for mutation and growth rates when fitting the multi-stage clonal expansion model to epidemiological IBD-CRC data yielded similar estimates to those found previously in sporadic CRC but suggest higher mutation rates and slightly lower growth rates in IBD. However, this data may not perfectly represent the natural history, as surveillance colonoscopy with lesion removal and colectomy alter the observable progression. Further, fitting to cancer incidence data alone presents parameter identifiability issues, restricting our initial fit to four parameters. To address these limitations, our study draws upon extensive clinical data from the U.S. Veterans Health Administration, employing validated methods using large language models to construct high-quality datasets with detailed information on surveillance colonoscopy timing, colectomies, and intermediate lesions extracted from pathology reports. To integrate these data, we developed a complementary fast simulation model, which will be released as an R package. This simulation model incorporates clinical interventions, such as colonoscopy with size-dependent lesion removal. Our combined analytical and simulation approach captures the complex precancerous evolution in IBD, providing a quantitative foundation for more effective, personalized surveillance guidelines. Further, this approach can be adapted to improve surveillance in the general population.

A biopsy is a common procedure used to diagnose diseases like cancer, infections, or inflammatory conditions. In cell population studies, biopsy samples provide valuable data to analyze cellular growth, proliferation rates, and structural abnormalities, which are essential for understanding disease progression. Estimating the growth (proliferation) rate of human cells is a challenging task. To address this, we have developed a method based on the birth-death Markov process to simulate the logistic growth model. We applied an extended Mean Field Approximation (MFA) for birth-death Markov processes, which accounts for fluctuations in the evolution of observables, such as moments. By calculating the theoretical moments from the birth-death process, we solved the inverse problem and estimated the growth rate. Additionally, we performed Markov Chain Monte Carlo (MCMC) simulations for both logistic growth and logistic growth with the Allee effect. The moments of the simulated population were used to predict the growth rate through regression analysis, achieving a high R-squared value. Finally, by applying this approach to biopsy data, one can estimate the proliferation rate of human cells with greater accuracy.

Introduction: Mathematical modeling is a key tool for understanding solid tumor biophysics, progression, and treatment resistance. Biophysical changes, such as elevated interstitial fluid pressure (IFP), are identified as major barriers to effective drug delivery. Incorporating patient-specific data into mathematical models offers the potential for personalized prognosis and treatment strategies for cancer patients. In this study, we explored the integration of patient-specific data from dynamic contrast-enhanced MRI (DCE-MRI) and diffusion-weighted MRI (DW-MRI) into a computational fluid dynamics (CFD) model of solid tumors to estimate the IFP and drug penetration profiles. Methods: As part of a translational study (EC/2019/1330, approved by Ghent University Hospital, Belgium), a patient with peritoneal metastasis underwent multi-sequential MRI, including T1-weighted (T1w) anatomical imaging, DCE-MRI, and DW-MRI. Tumor interstitial fluid pressure (IFP) was directly measured using a pressure transducer-tipped catheter for model validation. The CFD tumor model described interstitial fluid flow using Darcy’s law, the continuity equation, and Starling’s law, while drug penetration was modeled via the convection-diffusion-reaction equation. The 3D tumor geometry was derived from T1w images, vascular permeability from DCE-MRI using the extended-Tofts model, and hydraulic conductivity from DW-MRI. Results: An elevated IFP zone was observed in the central region of the tumor (up to 14 mmHg), while a lower IFP zone appeared at the tumor's edge. The clinically recorded IFP values (12.0 ± 2.5 mmHg) corresponded well with the simulation results. Drug penetration varied across the tumor, with deeper penetration in low-IFP regions. Conclusion: An image-based CFD model captured IFP and drug distribution variability, aligning with clinical data. This approach advances personalized oncology, potentially improving treatment strategies through noninvasive, patient-specific modeling.

The immune system can eradicate cancer, but various immunosuppressive mechanisms active within a tumor curb this beneficial response. Cytotoxic T cells (CD8s), regulatory T cells (Tregs), and antigen-presenting dendritic cells (DCs) play an important role in the immune response; however, it is very cumbersome to unravel the effects of multimodal interactions between tumor and immune cells and their contributions to tumor control using an experimental approach alone. Thus, to better understand the mechanisms that govern the interactions between immune cells and tumor cells and to identify the critical immunological features associated with tumor control and tumor escape, we built a mechanistic mathematical model of CD8s, Tregs, DCs, and tumor cells. The model accounts for tumor immunogenicity, the effects of IL-2 prolonging T cell lifespan, Treg suppression of antitumor immune response through CTLA-4, recruitment of immune cells into the tumor environment, and interferon-gamma upregulation of PD-L1 on DCs and tumor cells to deactivate T cells. We successfully fit the model to experimental data of tumor and immune cell dynamics. Employing Latin Hypercube Sampling, we generated over 1000 parameter sets that capture the sensitivity of αPD-1 immunotherapy. By comparing the parameter sets, we gain an insight into which mechanisms impinge the success of immunotherapy. We are now utilizing the model to explore combination immunotherapies that enhance the immune response in partial- and non-responders of αPD-1 immunotherapy. In particular, we are investigating what combinations of αPD-1, αCTLA-4, αICOSL, and αLAG-3 will stimulate CD8 activation without promoting Treg activation. After identifying the top combination therapies, we will validate our predictions experimentally. This integrated approach of modeling and experimental validation aims to advance our understanding of tumor-immune interactions and guide the development of more effective immunotherapeutic strategies.

Cancer is a significant global health issue, with treatment challenges arising from intratumor heterogeneity. This study examines the complex relationship between somatic evolution and phenotypic plasticity, explicitly focusing on the interplay between cell migration and proliferation [1]. We propose that evolution does not act directly on phenotypic traits, like the proliferation rate, but on the phenotypic plasticity in response to the microenvironment [2]. We study this hypothesis using a novel, spatially explicit model that tracks individual cells' phenotypic and genetic states. We assume cells change between mobile and growing states controlled by inherited and mutation-driven genotypes and the cells' microenvironment. We observe that cells at the tumor edge evolve to favor migration over proliferation and vice versa in the tumor bulk. However, this phenotypic heterogeneity can be realized by distinct regulations of the phenotypic switch, which depend on the apoptosis rate and the cells' ability to sense their environment. Emerging synthetic tumors display varying levels of heterogeneity, which we show are predictors of the cancer's recurrence time after treatment. Interestingly, higher phenotypic heterogeneity predicts poor treatment outcomes, unlike genetic heterogeneity. [1] Hatzikirou, H. et al. (2010). 'Go or Grow': the key to the emergence of invasion in tumour progression? Math. Med. Biol., 29(1), 49-65. [2] Syga S. et al. (2024) Evolution of phenotypic plasticity leads to tumor heterogeneity with implications for therapy. PLOS Comput. Biol. 20(8): e1012003.

A pivotal aspect of developing effective immunotherapies for solid tumors is the robust testing of product efficacy inside in vitro platforms.Collaborating with an experimental team that developed a novel microfluidic device at Children’s National Hospital (CNH), we developed a mathematical model to investigate immune cell migration and cytotoxicity within the device. Specifically, we study Chimeric Antigen Receptor (CAR) T-cell migration inside the channels, treating the cell as a moving boundary driven by a chemoattractant concentration gradient. The chemoattractant concentration is governed by two partial differential equations (PDEs) that incorporate key geometric elements of the device. We examine the motion of the cell as a function of its occlusion of the channel and find that certain cell shapes allow for multiple cells to travel inside the channel simultaneously. Additionally, we identify parameter regimes under which cells clog the channel, impairing their movement. All our findings are validated against experimental data provided by CNH. We integrate our model results into a broader model of the device, which also examines the cytotoxicity of CAR T-cells. This provides a tool for distinguishing experimental artifacts from genuine CAR T-cell behavior. This collaboration enabled the team at Children’s National Hospital to refine experimental conditions and uncover mechanisms enhancing CAR T-cell efficacy. [1] D Irimia, G Charras, N Agrawal, T Mitchison, M Toner, Polar stimulation and constrained cell migration in microfluidic channels,, Lab on a Chip 7 (12), 1783-1790

Extrachromosomal DNA (ecDNA) oncogene amplification is associated with treatment resistance and shorter survival in cancer. Currently, the spatial dynamics of ecDNA, and their evolutionary impact, are poorly understood. Here, we investigate ecDNA spatial-temporal evolution by integrating computational modeling with samples from 94 treatment-naive human IDH-wildtype glioblastoma patients. Random ecDNA segregation combined with ecDNA-conferred fitness advantages induce predictable spatial ecDNA copy-number patterns which depend on ecDNA oncogenic makeup. EGFR-ecDNAs often reach high copy-number, confer strong fitness advantages and do not co-amplify other oncogenes on the same ecDNA. In contrast, PDGFRA-ecDNAs reach lower copy-number, confer weaker fitness advantages and co-amplify other oncogenes. EGFR structural variants occur exclusively on ecDNA, arise from and are intermixed with wild-type EGFR-ecDNAs. Modeling suggests wild-type and variant EGFR-ecDNAs often accumulate before clonal expansion, even in patients co-amplifying multiple ecDNA species. Early emergence of oncogenic ecDNA under strong positive selection is confirmed in vivo and in vitro in mouse neural stem cells. Our results implicate ecDNA as a driver of gliomagenesis, and suggest a potential time window in which early ecDNA detection may facilitate more effective intervention.

Monoclonal antibody (mAb)-based therapeutics are pivotal in treating a wide range of diseases, including cancer. One key mechanism by which these antibodies exert anti-tumour effects is through antibody-dependent cellular cytotoxicity (ADCC). In ADCC, mAbs bind to specific antigens on tumour cells and Fc receptors on immune effector cells. This trimeric complex triggers these effector cells to kill the tumour. ADCC is influenced by multiple factors, notably the properties of the mAb and its interactions with Fc receptors and target antigens. However, the optimum conditions for ADCC remain unclear. In this study, we investigate how variations in target antigen and mAb properties, particularly antibody valency, modulate ADCC response to identify parameters that maximize its potency. We developed an ordinary differential equation (ODE) model to simulate mAb binding within the immune synapse and quantify trimeric complex formation. To link the number of trimeric complexes to ADCC response, we validated the model using Bayesian inference on ADCC assay data. The results suggest that lower-affinity mAbs enhance ADCC by increasing the number of target cell-bound antibodies. Our validated model indicates that a “steric penalty” is necessary for bivalently target-bound versus monovalently target-bound antibodies. Due to constraints from dual antigen binding, these antibodies experience limited mobility, reducing Fc receptor engagement. After model validation, we explored variations in target expression, binding affinity, and antibody valency on ADCC potency, quantified by EC50. Our key finding is that the optimal binding affinity for maximizing ADCC potency depends on antibody valency. Monovalent antibodies are most potent at high affinity, while bivalent antibodies peak at lower affinities. Furthermore, the magnitude of this effect varies with target expression levels.

Extrachromosomal DNA (ecDNA) oncogene amplification is associated with treatment resistance and shorter survival in cancer. Currently, the spatial dynamics of ecDNA, and their evolutionary impact, are poorly understood. Here, we investigate ecDNA spatial-temporal evolution by integrating computational modeling with samples from 94 treatment-naive human IDH-wildtype glioblastoma patients. Random ecDNA segregation combined with ecDNA-conferred fitness advantages induce predictable spatial ecDNA copy-number patterns which depend on ecDNA oncogenic makeup. EGFR-ecDNAs often reach high copy-number, confer strong fitness advantages and do not co-amplify other oncogenes on the same ecDNA. In contrast, PDGFRA-ecDNAs reach lower copy-number, confer weaker fitness advantages and co-amplify other oncogenes. EGFR structural variants occur exclusively on ecDNA, arise from and are intermixed with wild-type EGFR-ecDNAs. Modeling suggests wild-type and variant EGFR-ecDNAs often accumulate before clonal expansion, even in patients co-amplifying multiple ecDNA species. Early emergence of oncogenic ecDNA under strong positive selection is confirmed in vivo and in vitro in mouse neural stem cells. Our results implicate ecDNA as a driver of gliomagenesis, and suggest a potential time window in which early ecDNA detection may facilitate more effective intervention.

Monoclonal antibody (mAb)-based therapeutics are pivotal in treating a wide range of diseases, including cancer. One key mechanism by which these antibodies exert anti-tumour effects is through antibody-dependent cellular cytotoxicity (ADCC). In ADCC, mAbs bind to specific antigens on tumour cells and Fc receptors on immune effector cells. This trimeric complex triggers these effector cells to kill the tumour. ADCC is influenced by multiple factors, notably the properties of the mAb and its interactions with Fc receptors and target antigens. However, the optimum conditions for ADCC remain unclear. In this study, we investigate how variations in target antigen and mAb properties, particularly antibody valency, modulate ADCC response to identify parameters that maximize its potency. We developed an ordinary differential equation (ODE) model to simulate mAb binding within the immune synapse and quantify trimeric complex formation. To link the number of trimeric complexes to ADCC response, we validated the model using Bayesian inference on ADCC assay data. The results suggest that lower-affinity mAbs enhance ADCC by increasing the number of target cell-bound antibodies. Our validated model indicates that a “steric penalty” is necessary for bivalently target-bound versus monovalently target-bound antibodies. Due to constraints from dual antigen binding, these antibodies experience limited mobility, reducing Fc receptor engagement. After model validation, we explored variations in target expression, binding affinity, and antibody valency on ADCC potency, quantified by EC50. Our key finding is that the optimal binding affinity for maximizing ADCC potency depends on antibody valency. Monovalent antibodies are most potent at high affinity, while bivalent antibodies peak at lower affinities. Furthermore, the magnitude of this effect varies with target expression levels.

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

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

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

Musculoskeletal injuries can significantly impact the careers of racehorses, and are a common cause of lost training days and fatality. Most bone injuries arise from the gradual accumulation of microcracks through repeated training rather than spontaneous events. If severe enough, cracks may propagate through the bone and lead to fractures, often necessitating euthanasia. While training promotes bone adaptation to higher mechanical loads, overtraining can cause excessive damage. To better understand the biological processes underlying bone injury, we developed a lumped parameter model that combines the processes of bone adaptation, microdamage accumulation and bone repair. The model parameters were calibrated to experimental observations of bone volume fraction and time to fracture in racehorses. A sensitivity analysis identified joint loads (arising from training speed) and strides per day (arising from training distance) as key factors contributing to bone damage. Simulations of a typical training program showed that the majority of damage is incurred from training at racing speeds. While some microdamage is repaired during training, our model estimates that this process is insufficient to offset the damage accumulated. These findings emphasise the critical role of regular rest in preventing bone injury.

Diabetic kidney disease (DKD) progression is often marked by early glomerular endothelial cell (GEC) dysfunction, including alterations in the fenestration size and number linked to impaired glomerular filtration. However, the cellular mechanisms regulating GEC fenestrations remain poorly understood due to limitations in existing in vitro models, challenges in imaging small fenestrations in vivo, and inconsistencies between in vitro and in vivo findings. This study used a logic-based protein-protein interaction network model with normalized Hill functions for dynamics to explore how glucose-mediated signaling dysregulation impacts fenestration dynamics in GECs. We identified key drivers of fenestration loss and size changes by incorporating signaling pathways related to actin remodeling, myosin light chain kinase, Rho-associated kinase, calcium, and VEGF and its receptor. The model predicted how hyperglycemia in diabetic mice leads to significant fenestration loss and increased size of fenestrations. We found that glycemic control in the pre-DKD stage mitigated signaling dysregulation but was less effective as DKD developed and progressed. The model suggested alternative disease intervention strategies to maintain fenestration structure integrity, such as targeting Rho-associated kinase, VEGF-A, NFκB, and actin stress fibers.

Mutation accumulation (MA) experiments are crucial for understanding evolution. In microbial populations, these experiments typically involve periods of population growth, where a single individual forms a visible colony, followed by severe bottlenecks. Studies on the effects of positive and negative selection in MA experiments have shown that, for example, with 20 generations of growth between bottlenecks, beneficial mutations will be substantially over-represented cite{}; this effect is known as ``selection bias''. In previous work, we developed a fully stochastic discrete-time model that includes realistic offspring distributions, accounting for genetic drift and allowing for the loss of rare lineages. We demonstrated that when drift is included, selection bias is even stronger than previously predicted cite{}. Since bacterial division is unlikely to remain synchronized over 20 or more generations, this study extends the discrete-time model to a continuous-time framework. Since lineages that start reproducing early accrue a compounded advantage, a continuous-time model offers an even more accurate correction for selection in MA experiments.