Journal of Theoretical Biology

It has been hypothesised that the structure of tissues and the hierarchy of differentiation from stem cell to terminally-differentiated cell play a significant role in reducing the incidence of cancer in that tissue. One specific mechanism by which this risk can be reduced is by minimising the number of divisions - and hence the mutational risk - that cells accumulate as they divide to maintain tissue homeostasis. Here we investigate a mathematical model of cell division in a hierarchical tissue, calculating and minimising the divisional load while constraining parameters such that homeostasis is maintained. We show that the minimal divisional load is achieved by binary division tress with progenitor cells incapable of selfrenewal. Contrary to the protection hypothesis, we find that an increased stem cell turnover can lead to lower divisional load. Furthermore, we find that the optimal tissue structure depends on the time horizon of the duration of homeostasis, with faster stem cell division favoured in short-lived organisms and more progenitor compartments favoured in longer-lived organisms.

Atherosclerotic plaques are fatty deposits in arterial walls and a major cause of heart attacks and strokes. Macrophage proliferation triggers plaque growth and instability, but the exact conditions that cause stable plaques to become unstable remain unclear. To provide an insight into the conditions for this transition, we apply bifurcation analysis to the lipid-structured atherosclerosis model proposed by Chambers et al. (Bull Math Biol 86(8):104, 2024). Our main contribution is that the reduced dynamics of the system remain meaningful even beyond previously identified limits of validity. Furthermore, along with numerical bifurcation methods, the use of fast-slow analysis, combined with Fenichels theory, identifies a saddle-node bifurcation at infinity. A sharp threshold exists where macrophage proliferation balances emigration. Below this balance, the system stabilises in a biologically reasonable state; contrary to above it, macrophage numbers and lipid load grow unboundedly, triggering instability and runaway inflammation. Trends in determinant and eigenvalues also support this threshold. Parameter scans and heatmaps demonstrate that increased proliferation or reduced emigration enhances the number of macrophages and the lipid content of the necrotic core. Efferocytosis rate modulates downstream severity but does not shift the primary threshold. These findings reconcile conflicting results on macrophage proliferation, demonstrating that it is protective when emigration sufficiently balances this process. In other words, co-targeting reduced macrophage proliferation and enhanced emigration could help maintain plaque stability and reduce the risk of acute cardiovascular events. While this remains a theoretical recommendation, it offers a potential therapeutic strategy that authorises further investigation in experimental and clinical settings.

A new model based on a dynamical equation for the virus to microbe ratio (VMR) during log phase population growth shows that an optimal balance occurs between a short lysis time with low burst size vs. a long lysis time with large burst size. The model predicts that interventions lowering phage adsorption by killing free virus and/or limiting their access to bacteria favors the evolution of an increased lysis time and higher burst per infecting microbe until the intervention either drives a virulent phage extinct or, for temperate phage, drives the phage from its lytic phase into its lysogenic phase. The model also predicts that along an environmental gradient of increasing primary productivity the optimal lysis time shortens along the gradient, implying that the lytic life cycle goes around faster along the gradient. ImportanceA new approach to modeling phage life history predicts that virus respond to interventions that limit their adsorption onto bacteria by evolving a longer lysis time. The new model also predicts that lysis time of virus in nature shortens and the virus life cycle goes around faster as environmental conditions favoring virus production increase. These predictions show that virus life-history traits are not arbitrary and can be predicted in advance based on environmental conditions.

Mutations arising during cancer evolution are typically categorized as either drivers or passengers, depending on whether they increase the cell fitness. Recently, McFarland et al. introduced the Tug-of-War model for the joint effect of rare advantageous drivers and frequent but deleterious passengers. We examine this model under two common but distinct frameworks, the Moran model and the branching process. We show that frequently used statistics are similar between a version of the Moran model and the branching process conditioned on the final cell count, under different selection scenarios. We infer the selection coefficients for three breast cancer samples, resulting in good fits of the shape of their Site Frequency Spectra. All fitted values for the selective disadvantage of passenger mutations are nonzero, supporting the view that they exert deleterious selection during tumorigenesis that driver mutations must compensate.

The study of hematopoietic stem cell (HSCs) maintenance and differentiation to supply the hematopoietic system presents unique challenges, given the complex regulation of the process and the difficulty in observing cellular interactions in the stem cell niche. Quantitative methods and tools have emerged as valuable mechanisms to address this issue; however, the stochasticity of HSCs presents significant challenges for mathematical modeling, especially when bridging the gap between theoretical models and experimental validation. In this work, we have built a flexible and user-friendly stochastic dynamical and spatial model for long-term HSCs (LT-HSCs) and short-term HSCs (ST-HSCs) that captures experimentally observed cellular variability and heterogeneity. Our model implements the behavior of LT-HSCs and ST-HSCs and predicts their homeostatic dynamics. Furthermore, our model can be modified to explore various biological scenarios, such as stress-induced perturbations mediated by apoptosis, and successfully implement these conditions. Finally, the model incorporates spatial dynamics, simulating cell behavior in a 2D environment by combining Brownian motion with spatially graded parameters. *Summary StatementThis study addresses the challenge of characterizing hematopoietic stem cell (HSC) dynamics by developing a flexible, user-friendly stochastic spatial model of long-term and short-term HSCs. The model captures observed cellular variability and heterogeneity, predicts homeostatic dynamics, can be adapted to simulate stress-induced perturbations like apoptosis, and incorporates a spatial component to analyze HSC movement within a bone marrow niche.

The widely accepted model of tumor growth assumes tumors grow exponentially at a constant rate during early tumorigenesis when populations are small. The possibility that tumors might exhibit altered slower growth dynamics or even net cell death below a critical tumor size has yet to be fully explored. Deterministic growth models are capable of describing larger populations because population variation becomes small compared with the average, but when the population being modeled is small, the inherent stochasticity of the birth and death process produces significant variation. Recent advances in high throughput data collection allow for precise and sufficiently large data sets needed to capture this variation. Therefore, we present a stochastic modeling framework to describe and test the potential for altered growth dynamics at small tumor populations.

Temperate phage can initiate lysis or lysogeny after infecting a bacterial host. The genetic switch between lysis and lysogeny is mediated by phage regulatory genes as well as host and environmental factors. Recently, a new class of decision switches was identified in phage of the SPbeta group, mediated by the extracellular release of small, phage-encoded peptides termed arbitrium. Arbitrium peptides can be taken up by bacteria prior to infection, modulating the decision switch in the event of a subsequent phage infection. Increasing concentration of arbitrium increases the chance that a phage infection will lead to lysogeny, rather than lysis. Although prior work has centered on the molecular mechanisms of arbitrium-induced switching, here we focus on how selective pressures impact the benefits of plasticity in switching responses. In this work, we examine the possible advantages of near-term adaptation of communication-based decision switches used by the SPbeta-like group. We combine a nonlinear population model with a control theoretic approach to evaluate the relationship between a putative phage reaction norm (i.e., the probability of lysogeny as a function of arbitrium) and the near-term time horizon. We show the adaptive potential of communication-based lysis-lysogeny responses and find that optimal switching between lysis to lysogeny increases near-term fitness compared to fixed responses. We further find that plastic responses are robust to the inclusion of cellular-level stochasticity. These findings provide a principled basis to explore the long-term evolution of phage-encoded decision systems mediated by extracellular decision-signaling molecules, and the feedback between phage reaction norms and ecological context.

AO_SCPLOWBSTRACTC_SCPLOWFirst-line cancer treatment frequently fails due to initially rare therapeutic resistance. An important clinical question is then how to schedule subsequent treatments to maximize the probability of tumour eradication. Here, we provide a theoretical solution to this problem by using mathematical analysis and extensive stochastic simulations within the framework of evolutionary rescue theory to determine how best to exploit the vulnerability of small tumours to stochastic extinction. Whereas standard clinical practice is to wait for evidence of relapse, we confirm a recent hypothesis that the optimal time to switch to a second treatment is when the tumour is close to its minimum size before relapse, when it is likely undetectable. This optimum can lie slightly before or slightly after the nadir, depending on tumour parameters. Given that this exact time point may be difficult to determine in practice, we study windows of high extinction probability that lie around the optimal switching point, showing that switching after the relapse has begun is typically better than switching too early. We further reveal how treatment dose and tumour demographic and evolutionary parameters influence the predicted clinical outcome, and we determine how best to schedule drugs of unequal efficacy. Our work establishes a foundation for further experimental and clinical investigation of this evolutionarily-informed "extinction therapy" strategy.

In predicting metastatic potential and improving treatment outcomes in cancer research, it is crucial that we understand the dynamics of cancer cell dormancy and reactivation. In this paper we propose, study, and evaluate a cancer growth model that incorporates cell death, dormancy, reactivation, and proliferation in the secondary sites. Using experimental data on mice, we test various statistical distributions and identify models that represent the asymmetry and variability observed in dormancy durations. Notably, the estimated cancer cell death rate remained consistent across all tested distributions, supporting its biological relevance as a robust parameter for modelling dormancy survival dynamics. The most suitable among the distributions we studied, the most suitable ones exhibit heavy-tails and asymmetric skewness; this aligns with the prolonged and rare dormancy periods expected of cancer cells. Our findings stress the importance of selecting appropriate statistical models for dormancy, both in predicting cancer cell reactivation events, and in informing therapeutic strategies that focus dormancy-driven metastasis.

Tree populations worldwide are facing an unprecedented threat from a variety of tree diseases and invasive pests. Their spread, exacerbated by increasing globalisation and climate change, has an enormous environmental, economic and social impact. Computational agent-based models are a popular tool for describing and forecasting the spread of tree diseases due to their flexibility and ability to reveal collective behaviours. In this paper we present a versatile agentbased model with a Gaussian infectivity kernel to describe the spread of a generic tree disease through a synthetic treescape. We then explore several methods of calculating the basic reproduction number R0, a characteristic measurement of disease infectivity, defining the expected number of new infections resulting from one newly infected individual throughout their infectious period. It is a useful comparative summary parameter of a disease and can be used to explore the threshold dynamics of epidemics through mathematical models. We demonstrate several methods of estimating R0 through the agent-based model, including contact tracing, inferring the Kermack-McKendrick SIR model parameters using the linear noise approximation, and an analytical approximation. As an illustrative example, we then use the model and each of the methods to calculate estimates of R0 for the ash dieback epidemic in the UK.

Ageing is driven by damage accumulation leading to a decline in function over time. In single-cell systems, in addition to this damage accumulation within individuals, asymmetric damage partitioning at cell division can play a crucial role in shaping demographic ageing patterns. Despite empirical single-cell studies providing quantitative data at the molecular and demographic level, a comprehensive theory of how cellular damage production and asymmetric partitioning propagate and influence demographic patterns is still lacking. Here, we present a generic and flexible damage model using a stochastic differential equation approach that incorporates stochastic damage accumulation and asymmetric damage partitioning at cell divisions. We formulate an analytical approximation linking cellular and damage parameters to demographic ageing patterns. Interestingly, the lifespan of cells follows an inverse-gaussian distribution whose underlying properties derive from cellular and damage parameters. We demonstrate how stochasticity (noise) in damage production, asymmetry in damage partitioning, and division frequency shape lifespan distribution. Confronting the model to various empirical E.coli data reveals non-exponential scaling in mortality rates, a scaling that cannot be captured by classical Gompertz-Makeham models. Our findings provide a deep understanding of how fundamental processes contribute to cellular damage dynamics and generate demographic patterns. Our damage models generic nature offers a valuable framework for investigating ageing in diverse biological systems. SignificanceAsymmetries and randomness in cellular events play important roles in establishing the diversity at evolutionary and demographic scales. Looking at single-cells, ageing processes are influenced by stochastic damage accumulation and asymmetric damage partitioning at cell divisions. Utilising stochastic differential equations, we develop a cellular damage model that encapsulate both noisy damage accumulation within cells and asymmetric damage partitioning among cells when they divide. In doing so we bridge molecular stochastic processes, demographic fates of individuals, and population level demographic distributions. Our model aligns with empirical data from E.coli single-cell studies and advances our understanding compared to traditional demographic models. The generic and adaptable nature of our model paves the way for broader applications in ageing research across biological systems, highlighting the influences of stochasticity and asymmetry on cellular biology.

This research note investigates the optimal size of pools for pooled, COVID-19 testing when positive pools will be followed up by individual tests of pool members. Formulae for the optimum are derived and provided. The analysis indicates that O_LIoptimal pool sizes are unlikely to exceed about 20 individuals in realistic situations, C_LIO_LIoptimal pool size is influenced by prevalence in the population and the extent to which infection is clustered within the population, and, C_LIO_LIpools are most efficiently comprised of people with homogeneous risk, with heterogeneity across pools. C_LI

Cell size homeostasis in bacteria is a fundamental problem in systems biology, where cells maintain growth and division over many generations despite intrinsic fluctuations. Identifying the underlying control mechanism--whether division is triggered by reaching a critical size (sizer ) or by adding a fixed size increment (adder )--is essential for understanding this process. These two hypotheses are widely studied, yet there is no guarantee that either fully captures the true biological mechanism. More fundamentally, it has been unclear whether the control mechanism is statistically identifiable at all from lineage data. We address this question by developing a likelihood-based framework that explicitly accounts for threshold dynamics modeled as an Ornstein-Uhlenbeck process. Division timing is formulated as the first-passage-time (FPT) of this stochastic process to a time-dependent barrier. However, the FPT distribution lacks a closed-form analytical expression, preventing direct derivation of the maximum likelihood estimator (MLE). We overcome this challenge by training a neural network to approximate the FPT distribution and integrating it into the likelihood function, preserving analytical structure up to the FPT term. Simulations demonstrate that our method reliably distinguishes between sizer and adder mechanisms under realistic conditions where heuristic methods fail, providing the first evidence that the underlying control mechanism is identifiable. This hybrid analytical-machine learning approach provides a generalizable framework for studying stochastic threshold-based regulation in biological systems. Code reproducing the results is available at https://github.com/RonTeichner/newBacteria.

Prion diseases are lethal neurodegenerative disorders such as mad cow disease in bovines, chronic wasting disease in cervids, and Creutzfeldt-Jakob disease in humans. They are caused when the prion protein PrPC misfolds into PrPSc, which is capable of inducing further misfolding in healthy PrPC proteins. Recent in vivo experiments show that pharmacological chaperones can temporarily prevent this conversion by binding to PrPC molecules, and thus constitute a possible treatment. A second strategic approach uses interferons to decrease the concentration of PrPSc. In order to study the quantitative effects of these treatments on prion proliferation, we develop a model using a non-linear system of ordinary differential equations. By evaluating their efficacy and potency, we find that interferons act at lower doses and achieve greater prion decay rates. However, there are benefits in combining them with pharmacological chaperones in a two-fold therapy. This research is crucial to guide future prion experiments and inform potential treatment protocols.Competing Interest StatementThe authors have declared no competing interest.View Full Text

Hematopoietic Stem Cells (HSC) are the cells that give rise to all other blood cells and, as such, they are crucial in the healthy development of individuals. Wiskott-Aldrich Syndrome (WAS) is a severe disorder affecting the regulation of hematopoietic cells and is caused by mutations in the WASP gene. We consider data from a revolutionary gene therapy clinical trial, where HSC harvested from 3 WAS patients bone marrow have been edited and corrected using viral vectors. Upon re-infusion into the patient, the HSC multiply and differentiate into other cell types. The aim is to unravel the cell multiplication and cell differentiation process, which has until now remained elusive. This paper models the replenishment of blood lineages resulting from corrected HSC via a multivariate, density-dependent Markov process and develops an inferential procedure to estimate the dynamic parameters given a set of temporally sparsely observed trajectories. Starting from the master equation, we derive a system of non-linear differential equations for the evolution of the first- and second-order moments over time. We use these moment equations in a generalized method-of-moments framework to perform inference. The performance of our proposal has been evaluated by considering different sampling scenarios and measurement errors of various strengths using a simulation study. We also compared it to another state-of-the-art approach and found that our method is statistically more efficient. By applying our method to the Wiskott-Aldrich Syndrome gene therapy data we found strong evidence for a myeloid-based developmental pathway of hematopoietic cells where fates of lymphoid and myeloid cells remain coupled even after the loss of erythroid potential. All code used in this manuscript can be found in the online Supplement, and the latest version of the code is available at github. com/dp3ll1n/SLCDP_v1.0.

Thrombopoietin (TPO) is the ligand of the Mpl receptor and the key regulator of megakaryopoiesis and platelet production. A loss or gain of the TPO-receptor function affects haematopoiesis and results in severe diseases in humans. Appropriate mouse strains are available to mimic both myeloproliferative neoplasm (MPN) and congenital amegakaryocytic thrombocytopenia (CAMT) resulting from TPO overexpression or knockout of TPO receptor Mpl on megakaryocytes and platelets, respectively. However, at a quantitative level it is not understood, how the known regulations can establish the impaired but stable disease phenotypes. Starting out from an established mathematical model for megakaryopoiesis, we aim to adapt it to both the healthy situation and to distinct diseased phenotypes mimicking MPN. We thereby identify, that some of the model parameters are invariant with respect to the mouse strain while others have to be estimated in a strain-dependent manner. A systematic process of parameter identification provides strong evidence that the well-known excess production of megakaryocytes and early progenitors in MPN is either directly contingent on Mpl expression of platelets and megakaryocytes or, alternatively, that the knockout of Mpl is not as precisely restricted to megakaryocytes and platelets but may also effect their progenitors. In conclusion, our analysis hints towards an opaque control mechanism rendering megakaryopoiesis at a yet unknown level, and awaiting further experimental evaluation.

A phylogenetic tree has three types of attributes: size, shape (topology), and branch lengths. Phylody-namic studies are often motivated by questions regarding the size of clades, nevertheless, nearly all of the inference methods only make use of the other two attributes. In this paper, we ask whether there is additional information if we consider tree size more explicitly in phylodynamic inference methods. To address this question, we first needed to be able to compute the expected tree size distribution under a specified phylodynamic model; perhaps surprisingly, there is not a general method for doing so -- it is known what this is under a Yule or constant rate birth-death model but not for the more complicated scenarios researchers are often interested in. We present three different solutions to this problem: using i) the deterministic limit; ii) master equations; and iii) an ensemble moment approximation. Using simulations, we evaluate the accuracy of these three approaches under a variety of scenarios and alternative measures of tree size (i.e., sampling through time or only at the present; sampling ancestors or not). We then use the most accurate measures for the situation, to investigate the added informational content of tree size. We find that for two critical phylodynamic questions -- i) is diversification diversity dependent? and, ii) can we distinguish between alternative diversification scenarios? -- knowing the expected tree size distribution under the specified scenario provides insights that could not be gleaned from considering the expected shape and branch lengths alone. The contribution of this paper is both a novel set of methods for computing tree size distributions and a path forward for richer phylodynamic inference into the evolutionary and epidemiological processes that shape lineage trees.

Understanding and predicting extinction risk is a central challenge in population biology. Mathematical models incorporating Allee thresholds are commonly used to understand population dynamics and to assess extinction risks. Inaccurate predictions can have serious consequences for conservation management. In this simulation study, we develop a likelihood-based inference and prediction workflow to estimate parameters, including the Allee threshold and population diffusivity parameters, using noisy count data generated using a well-defined discrete model. Although parameters are identifiable according to commonly used criteria, the accuracy of resulting predictions depends strongly on the quantity, quality, collection time and spatial resolution of the data. Our workflow demonstrates that seemingly reliable parameter estimates can lead to inaccurate predictions, highlighting the need for careful consideration of data quality and quantity to guide extinction-risk modelling and prediction. Open source software is provided on GitHub to replicate and extend all results considered.

Understanding the incubation period of Legionnaires disease is vital for accurate source-term identification. Traditionally, researchers estimate the dose-dependent incubation period from human outbreak data, but this method suffers from the inability to estimate the exposure dose retrospectively for each case. This challenge limit the precision of incubation-period analysis using human case data. Existing within-host models, such as ordinary differential equation (ODE)-based and discrete-event stochastic approaches, estimate the dose-dependent incubation period of Legionnaires disease. However, discrete-event models, while useful, are so computationally costly that the within-host dynamics must be simplified to solely the Legionella and macrophage interactions. This simplification makes the computation feasible, but precludes cytokine interactions and adaptive immune response modelling. In this paper, we develop a new approach to model the within-host dynamics of Legion-naires disease that focuses on reducing computational cost while maintaining accuracy. Specifically, we propose a hybrid framework that integrates and improves upon existing ODE and discrete event within-host models of Legionnaires disease. By integrating the previously developed ODE and discrete-event stochastic models with stochastic differential equation (SDE) models, we create a unified system that adapts dynamically throughout the infection process. We quantify the points at which each model becomes the optimal tool for describing the infection, resulting in a flexible simulation of disease dynamics. Our hybrid model aligns with observed human incubation-period data and is the first framework of its kind in this context. This advancement offers a more robust platform for testing additional biological assumptions and improving our understanding of Legionnaires disease.

Influenza infection results from tightly coupled interactions between viral replication, host immune responses, and the emergence of clinical symptoms. While mathematical models have extensively characterized viral and immune dynamics, the mechanistic link between immune activity and disease severity remains poorly understood. Here, we develop an integrative within-host modeling framework that explicitly connects infection dynamics, immune responses, and symptom manifestation through a unified dynamical system. Using murine influenza data, we incorporate key immune components alongside a mechanistic representation of symptom progression, quantified via host weight loss. Our analysis identifies inflammatory signaling, particularly TNF--mediated pathways, as a central driver linking immune activity to symptom severity. Importantly, we demonstrate that age-dependent alterations in immune regulation reshape this coupling: aged hosts exhibit prolonged inflammatory responses that amplify and sustain symptom burden despite comparable viral kinetics. These results highlight that disease severity cannot be inferred from viral load alone, but instead emerges from the dynamical interplay between immune regulation and host physiology. This framework provides a quantitative basis for understanding age-specific morbidity and offers a foundation for designing interventions that target immune-mediated pathology rather than viral replication alone.