Approximate accelerated stochastic simulation of chemically reacting systems

Chiral symmetry breaking and information accumulation in pre-biological protocell evolution

We study a linear evolutionary model based on the two-dimensional distribution of protocells by total enantiomeric excess and the amount of stored information, which they can pass from generation to generation, and without any mutual inhibition. We show that the evolution of such systems occurs in four distinct stages. The first stage is an exponential growth of the concentration of protocells near the point [Formula: see text] and it should take negligible time on a geological scale. The second stage is a diffusion-like process in both dimensions. This process can also be accompanied by a drift in the direction of increased information passed from generation to generation, provided that the appropriate linear coefficient in the information storage subspace is large enough. The third stage is a rapid symmetry breaking and formation of the species near [Formula: see text] value of enantiomeric excess (assuming a small positive global enantiomeric asymmetry factor). The fourth stage is a relaxation toward a global stationary point, which is a narrow peak located near [Formula: see text] value of enantiomeric excess and some optimal value of the amount of stored information.

Efficient coupling of within-and between-host infectious disease dynamics

Mathematical models of infectious disease transmission typically neglect within-host dynamics. Yet within-host dynamics - including pathogen replication, host immune responses, and interactions with microbiota - are crucial not only for determining the progression of disease at the individual level, but also for driving within-host evolution and onwards transmission, and therefore shape dynamics at the population level. Various approaches have been proposed to model both within- and between-host dynamics, but these typically require considerable simplifying assumptions to couple processes at contrasting scales (e.g., the within-host dynamics quickly reach a steady state) or are computationally intensive. Here we propose a novel, readily adaptable and broadly applicable method for modelling both within- and between-host processes which can fully couple dynamics across scales and is both realistic and computationally efficient. By individually tracking the deterministic within-host dynamics of infected individuals, and stochastically coupling these to continuous host state variables at the population-level, we take advantage of fast numerical methods at both scales while still capturing individual transient within-host dynamics and stochasticity in transmission between hosts. Our approach closely agrees with full stochastic individual-based simulations and is especially useful when the within-host dynamics do not rapidly reach a steady state or over longer timescales to track pathogen evolution. By applying our method to different pathogen growth scenarios we show how common simplifying assumptions fundamentally change epidemiological and evolutionary dynamics.

Shared sanitation facilities and risk of respiratory virus transmission in resource-poor settings: A COVID-19 modeling case study

Water supply and sanitation are essential household services frequently shared in resource-poor settings. Shared sanitation can increase the risk of enteric pathogen transmission due to suboptimal cleanliness of facilities used by large numbers of individuals. It also can potentially increase the risk of respiratory disease transmission. As sanitation is an essential need, shared sanitation facilities may act as important respiratory pathogen transmission venues even with strict control measures such as stay-at-home recommendations in place. This analysis explores how behavioral and infrastructural conditions surrounding shared sanitation may individually and interactively influence respiratory pathogen transmission. We developed an individual-based community transmission model using COVID-19 as a motivating example parameterized from empirical literature to explore how transmission in shared latrines interacts with transmission at the community level. We explored mitigation strategies, including infrastructural and behavioral interventions. Our review of empirical literature confirms that shared sanitation venues in resource-poor settings are relatively small with poor ventilation and high use patterns. In these contexts, shared sanitation facilities may act as strong drivers of respiratory disease transmission, especially in areas reliant on shared facilities. Decreasing dependence on shared latrines was most effective at attenuating sanitation-associated transmission. Improvements to latrine ventilation and handwashing behavior were also able to decrease transmission. The type and order of interventions are important in successfully attenuating disease risk, with infrastructural and engineering controls being most effective when administered first, followed by behavioral controls after successful attenuation of sufficient alternate transmission routes. Beyond COVID-19, our modeling framework can be extended to address water, sanitation, and hygiene measures targeted at a range of environmentally mediated infectious diseases.

Tackling the Temporal Stiffness of Kinetic Monte Carlo Simulations of Well-Mixed Chemical Systems via On-the-Fly Scaling and Cost-Error Optimization

Reaction kinetics in biological systems are often subject to stochastic effects due to the low populations of reacting molecules, necessitating the adoption of kinetic Monte Carlo methods for their study. Such methods, however, can be computationally expensive, especially in the case of stiff systems, where some reactions are executed at much higher frequencies than others. We present an algorithm that reduces the reaction rate constants of the fast processes on-the-fly, thereby saving computational time, while keeping the approximation error within desirable limits. The algorithm couples the Modified Next Reaction Method for simulating stochastic systems with the Common Random Number framework and calculates accurate metrics for both the computational cost and approximation error by generating multiple sets of trajectories that correspond to increasingly reduced (downscaled) reaction rate constants. The optimum downscale factor is chosen via optimization of two conflicting objectives: (a) maximizing the speedup and (b) minimizing the approximation error introduced, and it is straightforward to tune the performance of the method, favoring accuracy versus speed or vice versa. Our approach is demonstrated on a biology-inspired well-mixed stiff system and is shown to accelerate the stochastic simulation thereof from 66 h down to 90 min, achieving a speed-up factor of 44×, without distorting the dynamics of the system studied.

Maintenance of long-term transposable element activity through regulation by nonautonomous elements

Transposable elements are DNA sequences that can move and replicate within genomes. Broadly, there are 2 types: autonomous elements, which encode the necessary enzymes for transposition, and nonautonomous elements, which rely on the enzymes produced by autonomous elements for their transposition. Nonautonomous elements have been proposed to regulate the numbers of transposable elements, which is a possible explanation for the persistence of transposition activity over long evolutionary times. However, previous modeling studies indicate that interactions between autonomous and nonautonomous elements usually result in the extinction of one type. Here, we study a stochastic model that allows for the stable coexistence of autonomous and nonautonomous retrotransposons. We determine the conditions for this coexistence and derive an analytical expression for the stationary distribution of their copy numbers, showing that nonautonomous elements regulate stochastic fluctuations and the number of autonomous elements in stationarity. We find that the stationary variances of each element can be expressed as a function of the average copy numbers and their covariance, enabling data comparison and model validation. These results suggest that continued transposition activity of transposable elements, regulated by nonautonomous elements, is a possible evolutionary outcome that could for example explain the long coevolutionary history of autonomous LINE1 and nonautonomous Alu element transposition in the human ancestry.

Selection for altruistic defense in structured populations

We model natural selection for or against an anti-parasite (or anti-predator) defense allele in a host (or prey) population that is structured into many demes. The defense behavior has a fitness cost for the actor compared to non defenders ("cheaters") in the same deme and locally reduces parasite growth rates. Hutzenthaler et al. (2022) have analytically derived a criterion for fixation or extinction of defenders in the limit of large populations, many demes, weak selection and slow migration. Here, we use both individual-based and diffusion-based simulation approaches to analyze related models. We find that the criterion still leads to accurate predictions for settings with finitely many demes and with various migration patterns. A key mechanism of providing a benefit of the defense trait is genetic drift due to randomness of reproduction and death events leading to between-deme differences in defense allele frequencies and host population sizes. We discuss an inclusive-fitness interpretation of this mechanism and present in-silico evidence that under these conditions a defense trait can be altruistic and still spread in a structured population.

Time-varying reproduction number estimation: fusing compartmental models with generalized additive models

The reproduction number, the mean number of secondary cases infected by each primary case, gives an indication of the effort required to control the disease. Beyond the well-known basic reproduction number, there are two natural extensions, namely the control and effective reproduction numbers. As behaviour, population immunity and viral characteristics can change with time, these reproduction numbers can vary over time. Real-world data can be complex, so in this work we consider a generalized additive model to smooth surveillance data through the explicit incorporation of day-of-the-week effects, to provide a simple measure of the time-varying growth rate associated with the data. Converting the resulting spline into an estimator for both the control and effective reproduction numbers requires assumptions on a model structure, which we here assume to be a compartmental model. The reproduction numbers calculated are based on both simulated and real-world data, and are compared with estimates from an already existing tool. The derived method for estimating the time-varying reproduction number is effective, efficient and comparable with other methods. It provides a useful alternative approach, which can be included as part of a toolbox of models, that is particularly apt at smoothing out day-of-the-week effects in surveillance.

A hybrid tau-leap for simulating chemical kinetics with applications to parameter estimation

We consider the problem of efficiently simulating stochastic models of chemical kinetics. The Gillespie stochastic simulation algorithm (SSA) is often used to simulate these models; however, in many scenarios of interest, the computational cost quickly becomes prohibitive. This is further exacerbated in the Bayesian inference context when estimating parameters of chemical models, as the intractability of the likelihood requires multiple simulations of the underlying system. To deal with issues of computational complexity in this paper, we propose a novel hybrid τ-leap algorithm for simulating well-mixed chemical systems. In particular, the algorithm uses τ-leap when appropriate (high population densities), and SSA when necessary (low population densities, when discrete effects become non-negligible). In the intermediate regime, a combination of the two methods, which uses the properties of the underlying Poisson formulation, is employed. As illustrated through a number of numerical experiments, the hybrid τ offers significant computational savings when compared with SSA without, however, sacrificing the overall accuracy. This feature is particularly welcomed in the Bayesian inference context, as it allows for parameter estimation of stochastic chemical kinetics at reduced computational cost.

Modelling plausible scenarios for the Omicron SARS-CoV-2 variant from early-stage surveillance

We used a spatially explicit agent-based model of SARS-CoV-2 transmission combined with spatially fine-grained COVID-19 observation data from Public Health Scotland to investigate the initial rise of the Omicron (BA.1) variant of concern. We evaluated plausible scenarios for transmission rate advantage and vaccine immune escape relative to the Delta variant based on the data that would have been available at that time. We also explored possible outcomes of different levels of imposed non-pharmaceutical intervention. The initial results of these scenarios were used to inform the Scottish Government in the early outbreak stages of the Omicron variant. Using the model with parameters fit over the Delta variant epidemic, some initial assumptions about Omicron transmission rate advantage and vaccine escape, and a simple growth rate fitting procedure, we were able to capture the initial outbreak dynamics for Omicron. We found that the modelled dynamics hold up to retrospective scrutiny. The modelled imposition of extra non-pharmaceutical interventions planned by the Scottish Government at the time would likely have little effect in light of the transmission rate advantage held by the Omicron variant and the fact that the planned interventions would have occurred too late in the outbreak's trajectory. Finally, we found that any assumptions made about the projected distribution of vaccines in the model population had little bearing on the outcome, in terms of outbreak size and timing. Instead, it was the landscape of prior immunity that was most important.

EpiFusion: Joint inference of the effective reproduction number by integrating phylodynamic and epidemiological modelling with particle filtering

Accurately estimating the effective reproduction number (Rt) of a circulating pathogen is a fundamental challenge in the study of infectious disease. The fields of epidemiology and pathogen phylodynamics both share this goal, but to date, methodologies and data employed by each remain largely distinct. Here we present EpiFusion: a joint approach that can be used to harness the complementary strengths of each field to improve estimation of outbreak dynamics for large and poorly sampled epidemics, such as arboviral or respiratory virus outbreaks, and validate it for retrospective analysis. We propose a model of Rt that estimates outbreak trajectories conditional upon both phylodynamic (time-scaled trees estimated from genetic sequences) and epidemiological (case incidence) data. We simulate stochastic outbreak trajectories that are weighted according to epidemiological and phylodynamic observation models and fit using particle Markov Chain Monte Carlo. To assess performance, we test EpiFusion on simulated outbreaks in which transmission and/or surveillance rapidly changes and find that using EpiFusion to combine epidemiological and phylodynamic data maintains accuracy and increases certainty in trajectory and Rt estimates, compared to when each data type is used alone. We benchmark EpiFusion's performance against existing methods to estimate Rt and demonstrate advances in speed and accuracy. Importantly, our approach scales efficiently with dataset size. Finally, we apply our model to estimate Rt during the 2014 Ebola outbreak in Sierra Leone. EpiFusion is designed to accommodate future extensions that will improve its utility, such as explicitly modelling population structure, accommodations for phylogenetic uncertainty, and the ability to weight the contributions of genomic or case incidence to the inference.

A modified variational approach to noisy cell signaling

Signaling in cells is full of noise and, hence, described with stochastic biochemical models. Thus, an efficient computation algorithm for these fluctuating reactions is much needed. Apart from the very popular Monte Carlo simulation, methods based on probability distributions are frequently desired due to their analytical tractability and possible numerical advantages in diverse circumstances, among which the variational approach is the most notable. In this paper, new basis functions are proposed to better depict possibly complex distribution profiles, and an extra regularization scheme is supplied to the variational equation to remove occasional degeneracy-induced singularities during the evolution. The new extension is applied to four typical biochemical reaction models and restores the Gillespie results accurately but with greatly reduced simulation time. This modified variational approach is expected to work in a wide range of cell signaling networks.

Using serosurveys to optimize surveillance for zoonotic pathogens

Zoonotic pathogens pose a significant risk to human health, with spillover into human populations contributing to chronic disease, sporadic epidemics, and occasional pandemics. Despite the widely recognized burden of zoonotic spillover, our ability to identify which animal populations serve as primary reservoirs for these pathogens remains incomplete. This challenge is compounded when prevalence reaches detectable levels only at specific times of year. In these cases, statistical models designed to predict the timing of peak prevalence could guide field sampling for active infections. Thus, we develop a general model that leverages routinely collected serosurveillance data to optimize sampling for elusive pathogens by predicting time windows of peak prevalence. Using simulated data sets, we show that our methodology reliably identifies times when pathogen prevalence is expected to peak. Then, we demonstrate an implementation of our method using publicly available data from two putative Ebolavirus reservoirs, straw-colored fruit bats (Eidolon helvum) and hammer-headed bats (Hypsignathus monstrosus). We envision our method being used to guide the planning of field sampling to maximize the probability of detecting active infections, and in cases when longitudinal data is available, our method can also yield predictions for the times of year that are most likely to produce future spillover events. The generality and simplicity of our methodology make it broadly applicable to a wide range of putative reservoir species where seasonal patterns of birth lead to predictable, but potentially short-lived, pulses of pathogen prevalence.

Computational Synthetic Biology Enabled through JAX: A Showcase

Mathematical modeling is indispensable in synthetic biology but remains underutilized. Tackling problems, from optimizing gene networks to simulating intracellular dynamics, can be facilitated by the ever-growing body of modeling approaches, be they mechanistic, stochastic, data-driven, or AI-enabled. Thanks to progress in the AI community, robust frameworks have emerged to enable researchers to access complex computational hardware and compilation. Previously, these frameworks focused solely on deep learning, but they have been developed to the point where running different forms of computation is relatively simple, as made possible, notably, by the JAX library. Running simulations at scale on GPUs speeds up research, which compounds enable larger-scale experiments and greater usability of code. As JAX remains underexplored in computational biology, we demonstrate its utility in three example projects ranging from synthetic biology to directed evolution, each with an accompanying demonstrative Jupyter notebook. We hope that these tutorials serve to democratize the flexible scaling, faster run-times, easy GPU portability, and mathematical enhancements (such as automatic differentiation) that JAX brings, all with only minor restructuring of code.

Tensor product algorithms for inference of contact network from epidemiological data

We consider a problem of inferring contact network from nodal states observed during an epidemiological process. In a black-box Bayesian optimisation framework this problem reduces to a discrete likelihood optimisation over the set of possible networks. The cardinality of this set grows combinatorially with the number of network nodes, which makes this optimisation computationally challenging. For each network, its likelihood is the probability for the observed data to appear during the evolution of the epidemiological process on this network. This probability can be very small, particularly if the network is significantly different from the ground truth network, from which the observed data actually appear. A commonly used stochastic simulation algorithm struggles to recover rare events and hence to estimate small probabilities and likelihoods. In this paper we replace the stochastic simulation with solving the chemical master equation for the probabilities of all network states. Since this equation also suffers from the curse of dimensionality, we apply tensor train approximations to overcome it and enable fast and accurate computations. Numerical simulations demonstrate efficient black-box Bayesian inference of the network.

Investigation of P. vivax elimination via mass drug administration: A simulation study

Plasmodium vivax is the most geographically widespread malaria parasite. P. vivax has the ability to remain dormant (as a hypnozoite) in the human liver and subsequently reactivate, which makes control efforts more difficult. Given the majority of P. vivax infections are due to hypnozoite reactivation, targeting the hypnozoite reservoir with a radical cure is crucial for achieving P. vivax elimination. Stochastic effects can strongly influence dynamics when disease prevalence is low or when the population size is small. Hence, it is important to account for this when modelling malaria elimination. We use a stochastic multiscale model of P. vivax transmission to study the impacts of multiple rounds of mass drug administration (MDA) with a radical cure, accounting for superinfection and hypnozoite dynamics. Our results indicate multiple rounds of MDA with a high-efficacy drug are needed to achieve a substantial probability of elimination. This work has the potential to help guide P. vivax elimination strategies by quantifying elimination probabilities for an MDA approach.

Adaptive tau-leaping methods for microscopic-lattice kinetic Monte Carlo simulations

Traditional Kinetic Monte Carlo (KMC) approaches, rooted in Gillespie's stochastic simulation algorithm, become computationally demanding in systems with a large range of timescales. The goal of this work is to propose and study new adaptive lattice-KMC time integration strategies for spatially non-uniform systems. To that end, two novel adaptive tau-leaping methods and their corresponding time integration strategies are developed based on the idea of the "n-fold" direct KMC method. These strategies allow for the simultaneous execution of multiple reactions, advancing time by adaptively selected coarse increments. We present numerical experiments comparing the proposed methods with existing approaches in a catalytic surface kinetics application involving ammonia decomposition.

Computation of random time-shift distributions for stochastic population models

Even in large systems, the effect of noise arising from when populations are initially small can persist to be measurable on the macroscale. A deterministic approximation to a stochastic model will fail to capture this effect, but it can be accurately approximated by including an additional random time-shift to the initial conditions. We present a efficient numerical method to compute this time-shift distribution for a large class of stochastic models. The method relies on differentiation of certain functional equations, which we show can be effectively automated by deriving rules for different types of model rates that arise commonly when mass-action mixing is assumed. Explicit computation of the time-shift distribution can be used to build a practical tool for the efficient generation of macroscopic trajectories of stochastic population models, without the need for costly stochastic simulations. Full code is provided to implement the calculations and we demonstrate the method on an epidemic model and a model of within-host viral dynamics.

Quantifying the risk of spillover reduction programs for human health

Reducing spillover of zoonotic pathogens is an appealing approach to preventing human disease and minimizing the risk of future epidemics and pandemics. Although the immediate human health benefit of reducing spillover is clear, over time, spillover reduction could lead to counterintuitive negative consequences for human health. Here, we use mathematical models and computer simulations to explore the conditions under which unanticipated consequences of spillover reduction can occur in systems where the severity of disease increases with age at infection. Our results demonstrate that, because the average age at infection increases as spillover is reduced, programs that reduce spillover can actually increase population-level disease burden if the clinical severity of infection increases sufficiently rapidly with age. If, however, immunity wanes over time and reinfection is possible, our results reveal that negative health impacts of spillover reduction become substantially less likely. When our model is parameterized using published data on Lassa virus in West Africa, it predicts that negative health outcomes are possible, but likely to be restricted to a small subset of populations where spillover is unusually intense. Together, our results suggest that adverse consequences of spillover reduction programs are unlikely but that the public health gains observed immediately after spillover reduction may fade over time as the age structure of immunity gradually re-equilibrates to a reduced force of infection.

Revisiting kinetic Monte Carlo algorithms for time-dependent processes: From open-loop control to feedback control

Simulating stochastic systems with feedback control is challenging due to the complex interplay between the system's dynamics and the feedback-dependent control protocols. We present a single-step-trajectory probability analysis to time-dependent stochastic systems. Based on this analysis, we revisit several time-dependent kinetic Monte Carlo (KMC) algorithms designed for systems under open-loop-control protocols. Our analysis provides a unified alternative proof to these algorithms, summarized into a pedagogical tutorial. Moreover, with the trajectory probability analysis, we present a novel feedback-controlled KMC algorithm that accurately captures the dynamics systems controlled by an external signal based on the measurements of the system's state. Our method correctly captures the system dynamics and avoids the artificial Zeno effect that arises from incorrectly applying the direct Gillespie algorithm to feedback-controlled systems. This work provides a unified perspective on existing open-loop-control KMC algorithms and also offers a powerful and accurate tool for simulating stochastic systems with feedback control.

Fitting Epidemic Models to Data: A Tutorial in Memory of Fred Brauer

Fred Brauer was an eminent mathematician who studied dynamical systems, especially differential equations. He made many contributions to mathematical epidemiology, a field that is strongly connected to data, but he always chose to avoid data analysis. Nevertheless, he recognized that fitting models to data is usually necessary when attempting to apply infectious disease transmission models to real public health problems. He was curious to know how one goes about fitting dynamical models to data, and why it can be hard. Initially in response to Fred's questions, we developed a user-friendly R package, fitode, that facilitates fitting ordinary differential equations to observed time series. Here, we use this package to provide a brief tutorial introduction to fitting compartmental epidemic models to a single observed time series. We assume that, like Fred, the reader is familiar with dynamical systems from a mathematical perspective, but has limited experience with statistical methodology or optimization techniques.

Computation Implemented by the Interaction of Chemical Reaction, Clustering, and De-Clustering of Molecules

A chemical reaction and its reaction environment are intrinsically linked, especially within the confines of narrow cellular spaces. Traditional models of chemical reactions often use differential equations with concentration as the primary variable, neglecting the density heterogeneity in the solution and the interaction between the reaction and its environment. We model the interaction between a chemical reaction and its environment within a geometrically confined space, such as inside a cell, by representing the environment through the size of molecular clusters. In the absence of fluctuations, the interplay between cluster size changes and the activation and inactivation of molecules induces oscillations. However, in unstable environments, the system reaches a fluctuating steady state. When an enzyme is introduced to this steady state, oscillations akin to action potential spike trains emerge. We examine the behavior of these spike trains and demonstrate that they can be used to implement logic gates. We discuss the oscillations and computations that arise from the interaction between a chemical reaction and its environment, exploring their potential for contributing to chemical intelligence.

Final nanoparticle size distribution under unusual parameter regimes

We explore the large-scale behavior of a stochastic model for nanoparticle growth in an unusual parameter regime. This model encompasses two types of reactions: nucleation, where n monomers aggregate to form a nanoparticle, and growth, where a nanoparticle increases its size by consuming a monomer. Reverse reactions are disregarded. We delve into a previously unexplored parameter regime. Specifically, we consider a scenario where the growth rate of the first newly formed particle is of the same order of magnitude as the nucleation rate, in contrast to the classical scenario where, in the initial stage, nucleation dominates over growth. In this regime, we investigate the final size distribution as the initial number of monomers tends to infinity through extensive simulation studies utilizing state-of-the-art stochastic simulation methods with an efficient implementation and supported by high-performance computing infrastructure. We observe the emergence of a deterministic limit for the particle's final size density. To scale up the initial number of monomers to approximate the magnitudes encountered in real experiments, we introduce a novel approximation process aimed at faster simulation. Remarkably, this approximating process yields a final size distribution that becomes very close to that of the original process when the available monomers approach infinity. Simulations of the approximating process further support the conjecture of the emergence of a deterministic limit.

Poissonian Cellular Potts Models Reveal Nonequilibrium Kinetics of Cell Sorting

Cellular Potts models are broadly applied across developmental biology and cancer research. We overcome limitations of the traditional approach, which reinterprets a modified Metropolis sampling as ad hoc dynamics, by introducing a physical timescale through Poissonian kinetics and by applying principles of stochastic thermodynamics to separate thermal and relaxation effects from athermal noise and nonconservative forces. Our method accurately describes cell-sorting dynamics in mouse-embryo development and identifies the distinct contributions of nonequilibrium processes, e.g., cell growth and active fluctuations.

Advanced methods for gene network identification and noise decomposition from single-cell data

Central to analyzing noisy gene expression systems is solving the Chemical Master Equation (CME), which characterizes the probability evolution of the reacting species' copy numbers. Solving CMEs for high-dimensional systems suffers from the curse of dimensionality. Here, we propose a computational method for improved scalability through a divide-and-conquer strategy that optimally decomposes the whole system into a leader system and several conditionally independent follower subsystems. The CME is solved by combining Monte Carlo estimation for the leader system with stochastic filtering procedures for the follower subsystems. We demonstrate this method with high-dimensional numerical examples and apply it to identify a yeast transcription system at the single-cell resolution, leveraging mRNA time-course experimental data. The identification results enable an accurate examination of the heterogeneity in rate parameters among isogenic cells. To validate this result, we develop a noise decomposition technique exploiting time-course data but requiring no supplementary components, e.g., dual-reporters.

An Agent-Based Model Shows How Mixed Infections Drive Multiyear Pathotype Dynamics in a Plant-Virus System

Mathematical models are widely used to understand the evolution and epidemiology of plant pathogens under a variety of scenarios. Here, we used this approach to analyze the effects of different traits intrinsic and extrinsic to plant-virus interactions on the dynamics of virus pathotypes in genetically heterogeneous plant-virus systems. For this, we propose an agent-based epidemiological model that includes epidemiologically significant pathogen life-history traits related to virulence, transmission, and survival in the environment and allows for integrating long- and short-distance transmission, primary and secondary infections, and within-host pathogen competition in mixed infections. The study focuses on the tobamovirus-pepper pathosystem. Model simulations allowed us to integrate pleiotropic effects of resistance-breaking mutations on different virus life-history traits into the net costs of resistance breaking, allowing for predictions on multiyear pathotype dynamics. We also explored the effects of two control measures, the use of host resistance and roguing of symptomatic plants, that modify epidemiological attributes of the pathogens to understand how their populations will respond to evolutionary pressures. One major conclusion points to the importance of pathogen competition within mixed-infected hosts as a component of the overall fitness of each pathogen that, thus, drives their multiyear dynamics.

Self-organization of PIP3 waves is controlled by the topology and curvature of cell membranes

Phosphatidylinositol (3,4,5)-trisphosphate (PIP3) is a signaling lipid on the plasma membrane that plays a fundamental role in cell signaling with a strong impact on cell physiology and diseases. It is responsible for the protruding edge formation, cell polarization, macropinocytosis, and other membrane remodeling dynamics in cells. It has been shown that the membrane confinement and curvature affects the wave formation of PIP3 and F-actin. But, even in the absence of F-actin, a complex self-organization of the spatiotemporal PIP3 waves is observed. In recent findings, we have shown that these waves can be guided and pinned on strongly bended Dictyostelium membranes caused by molecular crowding and curvature-limited diffusion. Based on these experimental findings, we investigate the spatiotemporal PIP3 wave dynamics on realistic three-dimensional cell-like membranes to explore the effect of curvature-limited diffusion, as observed experimentally. We use an established stochastic reaction-diffusion model with enzymatic Michaelis-Menten-type reactions that mimics the dynamics of Dictyostelium cells. As these cells mimic the three-dimensional shape and size observed experimentally, we found that the PIP3 wave directionality can be explained by a Hopf-like and a reverse periodic-doubling bifurcation for uniform diffusion and curvature-limited diffusion properties. Finally, we compare the results with recent experimental findings and discuss the discrepancy between the biological and numerical results.

A seven-step guide to spatial, agent-based modelling of tumour evolution

Spatial agent-based models are frequently used to investigate the evolution of solid tumours subject to localized cell-cell interactions and microenvironmental heterogeneity. As spatial genomic, transcriptomic and proteomic technologies gain traction, spatial computational models are predicted to become ever more necessary for making sense of complex clinical and experimental data sets, for predicting clinical outcomes, and for optimizing treatment strategies. Here we present a non-technical step by step guide to developing such a model from first principles. Stressing the importance of tailoring the model structure to that of the biological system, we describe methods of increasing complexity, from the basic Eden growth model up to off-lattice simulations with diffusible factors. We examine choices that unavoidably arise in model design, such as implementation, parameterization, visualization and reproducibility. Each topic is illustrated with examples drawn from recent research studies and state of the art modelling platforms. We emphasize the benefits of simpler models that aim to match the complexity of the phenomena of interest, rather than that of the entire biological system. Our guide is aimed at both aspiring modellers and other biologists and oncologists who wish to understand the assumptions and limitations of the models on which major cancer studies now so often depend.

Minimal Mechanisms of Microtubule Length Regulation in Living Cells

The microtubule cytoskeleton is responsible for sustained, long-range intracellular transport of mRNAs, proteins, and organelles in neurons. Neuronal microtubules must be stable enough to ensure reliable transport, but they also undergo dynamic instability, as their plus and minus ends continuously switch between growth and shrinking. This process allows for continuous rebuilding of the cytoskeleton and for flexibility in injury settings. Motivated by in vivo experimental data on microtubule behavior in Drosophila neurons, we propose a mathematical model of dendritic microtubule dynamics, with a focus on understanding microtubule length, velocity, and state-duration distributions. We find that limitations on microtubule growth phases are needed for realistic dynamics, but the type of limiting mechanism leads to qualitatively different responses to plausible experimental perturbations. We therefore propose and investigate two minimally-complex length-limiting factors: limitation due to resource (tubulin) constraints and limitation due to catastrophe of large-length microtubules. We combine simulations of a detailed stochastic model with steady-state analysis of a mean-field ordinary differential equations model to map out qualitatively distinct parameter regimes. This provides a basis for predicting changes in microtubule dynamics, tubulin allocation, and the turnover rate of tubulin within microtubules in different experimental environments.

Mathematical Modelling of Parasite Dynamics: A Stochastic Simulation-Based Approach and Parameter Estimation via Modified Sequential-Type Approximate Bayesian Computation

The development of mathematical models for studying newly emerging and re-emerging infectious diseases has gained momentum due to global events. The gyrodactylid-fish system, like many host-parasite systems, serves as a valuable resource for ecological, evolutionary, and epidemiological investigations owing to its ease of experimental manipulation and long-term monitoring. Although this system has an existing individual-based model, it falls short in capturing information about species-specific microhabitat preferences and other biological details for different Gyrodactylus strains across diverse fish populations. This current study introduces a new individual-based stochastic simulation model that uses a hybrid τ -leaping algorithm to incorporate this essential data, enhancing our understanding of the complexity of the gyrodactylid-fish system. We compare the infection dynamics of three gyrodactylid strains across three host populations. A modified sequential-type approximate Bayesian computation (ABC) method, based on sequential Monte Carlo and sequential importance sampling, is developed. Additionally, we establish two penalised local-linear regression methods (based on L1 and L2 regularisations) for ABC post-processing analysis to fit our model using existing empirical data. With the support of experimental data and the fitted mathematical model, we address open biological questions for the first time and propose directions for future studies on the gyrodactylid-fish system. The adaptability of the mathematical model extends beyond the gyrodactylid-fish system to other host-parasite systems. Furthermore, the modified ABC methodologies provide efficient calibration for other multi-parameter models characterised by a large set of correlated or independent summary statistics.

The impact of inaccurate assumptions about antibody test accuracy on the parametrisation and results of infectious disease models of epidemics

The parametrisation of infectious disease models is often done based on epidemiological studies that use diagnostic and serology tests to establish disease prevalence or seroprevalence in the population being modelled. During outbreaks of an emerging infectious disease, tests are often used, both for disease control and epidemiological studies, before studies evaluating their accuracy in the population have concluded, with assumptions made about accuracy parameters like sensitivity and specificity. In this simulation study, we simulated such an outbreak, based on the case study of COVID-19, and found that inaccurate parametrisation of infectious disease models due to assumptions about antibody test accuracy in a seroprevalence study can cause modelling results that inform public health decisions to be inaccurate; for example, in our simulation setup, assuming that antibody test specificity was 0.99 instead of 0.90 when it was in fact 0.90 led to an average relative difference of 0.78 in model-projected peak hospitalisations, even when test sensitivity and all other parameters were accurately characterised. We therefore suggest that methods to speed up test evaluation studies are vitally important in the public health response to an emerging outbreak.

Mild HIV-specific selective forces overlaying natural CD4+ T cell dynamics explain the clonality and decay dynamics of HIV reservoir cells

The latent reservoir of HIV persists for decades in people living with HIV (PWH) on antiretroviral therapy (ART). To determine if persistence arises from the natural dynamics of memory CD4+ T cells harboring HIV, we compared the clonal dynamics of HIV proviruses to that of memory CD4+ T cell receptors (TCRβ) from the same PWH and from HIV-seronegative people. We show that clonal dominance of HIV proviruses and antigen-specific CD4+ T cells are similar but that the field's understanding of the persistence of the less clonally dominant reservoir is significantly limited by undersampling. We demonstrate that increasing reservoir clonality over time and differential decay of intact and defective proviruses cannot be explained by mCD4+ T cell kinetics alone. Finally, we develop a stochastic model of TCRβ and proviruses that recapitulates experimental observations and suggests that HIV-specific negative selection mediates approximately 6% of intact and 2% of defective proviral clearance. Thus, HIV persistence is mostly, but not entirely, driven by natural mCD4+ T cell kinetics.

Competition between Mononuclear and Binuclear Copper Sites across Different Zeolite Topologies

A key challenge for metal-exchanged zeolites is the determination of metal cation speciation and nuclearity under synthesis and reaction conditions. Copper-exchanged zeolites, which are widely used in automotive emissions control and potential catalysts for partial methane oxidation, have in particular evidenced a wide variety of Cu structures that are observed to change with exposure conditions, zeolite composition, and topology. Here, we develop predictive models for Cu cation speciation and nuclearity in CHA, MOR, BEA, AFX, and FER zeolite topologies using interatomic potentials, quantum chemical calculations, and Monte Carlo simulations to interrogate this vast configurational and compositional space. Model predictions are used to rationalize experimentally observed differences between Cu-zeolites in a wide-body of literature, including nuclearity populations, structural variations, and methanol per Cu yields. Our results show that both topological features and commonly observed Al-siting biases in MOR zeolites increase the population of binuclear Cu sites, explaining the small population of mononuclear Cu sites observed in these materials relative to other zeolites such as CHA and BEA. Finally, we used a machine learning classification model to determine the preference to form mononuclear or binuclear Cu sites at different Al configurations in 200 zeolites in the international zeolite database. Model results reveal several zeolite topologies at extreme ends of the mononuclear vs binuclear spectrum, highlighting synthetic options for realization of zeolites with strong Cu nuclearity preferences.

Protocol for biomodel engineering of unilevel to multilevel biological models using colored Petri nets

Biological systems inherently span multiple levels, which can pose challenges in spatial representation for modelers. We present a protocol that utilizes colored Petri nets to construct and analyze biological models of systems, encompassing both unilevel and multilevel scenarios. We detail a modeling workflow exploiting the PetriNuts platform comprising a set of tools linked together via common file formats. We describe steps for modeling preparation, component-level modeling and analysis, followed by system-level modeling and analysis, and model use.

Identification of dynamical changes of rabies transmission under quarantine: Community-based measures towards rabies elimination

Quarantine has been long used as a public health response to emerging infectious diseases, particularly at the onset of an epidemic when the infected proportion of a population remains identifiable and logistically tractable. In theory, the same logic should apply to low-incidence infections; however, the application and impact of quarantine in low prevalence settings appears less common and lacks a formal analysis. Here, we present a quantitative framework using a series of progressively more biologically realistic models of canine rabies in domestic dogs and from dogs to humans, a suitable example system to characterize dynamical changes under varying levels of dog quarantine. We explicitly incorporate health-seeking behaviour data to inform the modelling of contact-tracing and exclusion of rabies suspect and probable dogs that can be identified through bite-histories of patients presenting at anti-rabies clinics. We find that a temporary quarantine of rabies suspect and probable dogs provides a powerful tool to curtail rabies transmission, especially in settings where optimal vaccination coverage is yet to be achieved, providing a critical stopgap to reduce the number of human and animal deaths due to rabid bites. We conclude that whilst comprehensive measures including sensitive surveillance and large-scale vaccination of dogs will be required to achieve disease elimination and sustained freedom given the persistent risk of rabies re-introductions, quarantine offers a low-cost community driven solution to intersectoral health burden.

Identifying Knowledge Gaps through the Systematic Review of Temperature-Driven Variability in the Competence of Aedes aegypti and Ae. albopictus for Chikungunya Virus

Temperature is a well-known effector of several transmission factors of mosquito-borne viruses, including within mosquito dynamics. These dynamics are often characterized by vector competence and the extrinsic incubation period (EIP). Vector competence is the intrinsic ability of a mosquito population to become infected with and transmit a virus, while EIP is the time it takes for the virus to reach the salivary glands and be expectorated following an infectious bloodmeal. Temperatures outside the optimal range act on life traits, decreasing transmission potential, while increasing temperature within the optimal range correlates to increasing vector competence and a decreased EIP. These relatively well-studied effects of other Aedes borne viruses (dengue and Zika) are used to make predictions about transmission efficiency, including the challenges presented by urban heat islands and climate change. However, the knowledge of temperature and chikungunya (CHIKV) dynamics within its two primary vectors-Ae. aegypti and Ae. albopictus-remains less characterized, even though CHIKV remains a virus of public-health importance. Here, we review the literature and summarize the state of the literature on CHIKV and temperature dependence of vector competence and EIP and use these data to demonstrate how the remaining knowledge gap might confound the ability to adequately predict and, thus, prepare for future outbreaks.

Error-Controlled Coarse-Graining Dynamics with Mean-Field Randomization

In order to comprehend the stochastic behavior of biological systems, it is essential to accurately infer the dynamics of chemical reaction networks. However, computation of the likelihood remains a bottleneck. In this study, we propose the mean-field randomization procedure as a means of efficiently generating error-controlled coarse-graining dynamics. The error is measured by mutual information between the generated trajectories and the coarse-graining procedure. We demonstrate that the exact dynamics can be recovered by resampling, which eliminates the correlation between the dynamics and the procedure. We developed three algorithms to efficiently generate exact or coarse-graining trajectories within a specified error range. By subjecting our algorithms to testing on chemical reaction systems of varying complexities and scales, we observe that they outperform existing state-of-the-art algorithms, and the efficiency of coarse-graining trajectory generation is only weakly dependent on system scales.

Inferring delays in partially observed gene regulation processes

MOTIVATION

Cell function is regulated by gene regulatory networks (GRNs) defined by protein-mediated interaction between constituent genes. Despite advances in experimental techniques, we can still measure only a fraction of the processes that govern GRN dynamics. To infer the properties of GRNs using partial observation, unobserved sequential processes can be replaced with distributed time delays, yielding non-Markovian models. Inference methods based on the resulting model suffer from the curse of dimensionality.

RESULTS

We develop a simulation-based Bayesian MCMC method employing an approximate likelihood for the efficient and accurate inference of GRN parameters when only some of their products are observed. We illustrate our approach using a two-step activation model: an activation signal leads to the accumulation of an unobserved regulatory protein, which triggers the expression of observed fluorescent proteins. With prior information about observed fluorescent protein synthesis, our method successfully infers the dynamics of the unobserved regulatory protein. We can estimate the delay and kinetic parameters characterizing target regulation including transcription, translation, and target searching of an unobserved protein from experimental measurements of the products of its target gene. Our method is scalable and can be used to analyze non-Markovian models with hidden components.

AVAILABILITY AND IMPLEMENTATION

Our code is implemented in R and is freely available with a simple example data at https://github.com/Mathbiomed/SimMCMC.

Global quantitative understanding of non-equilibrium cell fate decision-making in response to pheromone

Cell-cycle arrest and polarized growth are commonly used to characterize the response of yeast to pheromone. However, the quantitative decision-making processes underlying time-dependent changes in cell fate remain unclear. In this study, we conducted single-cell level experiments to observe multidimensional responses, uncovering diverse fates of yeast cells. Multiple states are revealed, along with the kinetic switching rates and pathways among them, giving rise to a quantitative landscape of mating response. To quantify the experimentally observed cell fates, we developed a theoretical framework based on non-equilibrium landscape and flux theory. Additionally, we performed stochastic simulations of biochemical reactions to elucidate signal transduction and cell growth. Notably, our experimental findings have provided the first global quantitative evidence of the real-time synchronization between intracellular signaling, physiological growth, and morphological functions. These results validate the proposed underlying mechanism governing the emergence of multiple cell fate states. This study introduces an emerging mechanistic approach to understand non-equilibrium cell fate decision-making in response to pheromone.

The effect of random virus failure following cell entry on infection outcome and the success of antiviral therapy

A virus infection can be initiated with very few or even a single infectious virion, and as such can become extinct, i.e. stochastically fail to take hold or spread significantly. There are many ways that a fully competent infectious virion, having successfully entered a cell, can fail to cause a productive infection, i.e. one that yields infectious virus progeny. Though many stochastic models (SMs) have been developed and used to estimate a virus infection's establishment probability, these typically neglect infection failure post virus entry. The SM presented herein introduces parameter [Formula: see text] which corresponds to the probability that a virion's entry into a cell will result in a productive cell infection. We derive an expression for the likelihood of infection establishment in this new SM, and find that prophylactic therapy with an antiviral reducing [Formula: see text] is at least as good or better at decreasing the establishment probability, compared to antivirals reducing the rates of virus production or virus entry into cells, irrespective of the SM parameters. We investigate the difference in the fraction of cells consumed by so-called extinct versus established virus infections, and find that this distinction becomes biologically meaningless as the probability of establishment approaches zero. We explain why the release of virions continuously over an infectious cell's lifespan, rather than as a single burst at the end of the cell's lifespan, does not result in an increased risk of infection extinction. We show, instead, that the number of virus released, not the timing of the release, affects infection establishment and associated critical antiviral efficacy.

Catalyst: Fast and flexible modeling of reaction networks

We introduce Catalyst.jl, a flexible and feature-filled Julia library for modeling and high-performance simulation of chemical reaction networks (CRNs). Catalyst supports simulating stochastic chemical kinetics (jump process), chemical Langevin equation (stochastic differential equation), and reaction rate equation (ordinary differential equation) representations for CRNs. Through comprehensive benchmarks, we demonstrate that Catalyst simulation runtimes are often one to two orders of magnitude faster than other popular tools. More broadly, Catalyst acts as both a domain-specific language and an intermediate representation for symbolically encoding CRN models as Julia-native objects. This enables a pipeline of symbolically specifying, analyzing, and modifying CRNs; converting Catalyst models to symbolic representations of concrete mathematical models; and generating compiled code for numerical solvers. Leveraging ModelingToolkit.jl and Symbolics.jl, Catalyst models can be analyzed, simplified, and compiled into optimized representations for use in numerical solvers. Finally, we demonstrate Catalyst's broad extensibility and composability by highlighting how it can compose with a variety of Julia libraries, and how existing open-source biological modeling projects have extended its intermediate representation.

Evolutionary safety of lethal mutagenesis driven by antiviral treatment

Nucleoside analogs are a major class of antiviral drugs. Some act by increasing the viral mutation rate causing lethal mutagenesis of the virus. Their mutagenic capacity, however, may lead to an evolutionary safety concern. We define evolutionary safety as a probabilistic assurance that the treatment will not generate an increased number of mutants. We develop a mathematical framework to estimate the total mutant load produced with and without mutagenic treatment. We predict rates of appearance of such virus mutants as a function of the timing of treatment and the immune competence of patients, employing realistic assumptions about the vulnerability of the viral genome and its potential to generate viable mutants. We focus on the case study of Molnupiravir, which is an FDA-approved treatment against Coronavirus Disease-2019 (COVID-19). We estimate that Molnupiravir is narrowly evolutionarily safe, subject to the current estimate of parameters. Evolutionary safety can be improved by restricting treatment with this drug to individuals with a low immunological clearance rate and, in future, by designing treatments that lead to a greater increase in mutation rate. We report a simple mathematical rule to determine the fold increase in mutation rate required to obtain evolutionary safety that is also applicable to other pathogen-treatment combinations.

Design and Simulation of a Multilayer Chemical Neural Network That Learns via Backpropagation

The design and implementation of adaptive chemical reaction networks, capable of adjusting their behavior over time in response to experience, is a key goal for the fields of molecular computing and DNA nanotechnology. Mainstream machine learning research offers powerful tools for implementing learning behavior that could one day be realized in a wet chemistry system. Here we develop an abstract chemical reaction network model that implements the backpropagation learning algorithm for a feedforward neural network whose nodes employ the nonlinear "leaky rectified linear unit" transfer function. Our network directly implements the mathematics behind this well-studied learning algorithm, and we demonstrate its capabilities by training the system to learn a linearly inseparable decision surface, specifically, the XOR logic function. We show that this simulation quantitatively follows the definition of the underlying algorithm. To implement this system, we also report ProBioSim, a simulator that enables arbitrary training protocols for simulated chemical reaction networks to be straightforwardly defined using constructs from the host programming language. This work thus provides new insight into the capabilities of learning chemical reaction networks and also develops new computational tools to simulate their behavior, which could be applied in the design and implementations of adaptive artificial life.

Approximate simulation of cortical microtubule models using dynamical graph grammars

Dynamical graph grammars (DGGs) are capable of modeling and simulating the dynamics of the cortical microtubule array (CMA) in plant cells by using an exact simulation algorithm derived from a master equation; however, the exact method is slow for large systems. We present preliminary work on an approximate simulation algorithm that is compatible with the DGG formalism. The approximate simulation algorithm uses a spatial decomposition of the domain at the level of the system's time-evolution operator, to gain efficiency at the cost of some reactions firing out of order, which may introduce errors. The decomposition is more coarsely partitioned by effective dimension (d= 0 to 2 or 0 to 3), to promote exact parallelism between different subdomains within a dimension, where most computing will happen, and to confine errors to the interactions between adjacent subdomains of different effective dimensions. To demonstrate these principles we implement a prototype simulator, and run three simple experiments using a DGG for testing the viability of simulating the CMA. We find evidence indicating the initial formulation of the approximate algorithm is substantially faster than the exact algorithm, and one experiment leads to network formation in the long-time behavior, whereas another leads to a long-time behavior of local alignment.

A Stochastic Mobility-Driven Spatially Explicit SEIQRD covid-19 Model with VOCs, Seasonality, and Vaccines

In this work, we extend our previously developed compartmental SEIQRD model for sars-cov-2 in Belgium. We introduce sars-cov-2 variants of concern, vaccines, and seasonality in our model, as their addition has proven necessary for modelling sars-cov-2 transmission dynamics during the 2020-2021 covid-19 pandemic in Belgium. The model is geographically stratified into eleven spatial patches (provinces), and a telecommunication dataset provided by Belgium's biggest operator is used to incorporate interprovincial mobility. We calibrate the model using the daily number of hospitalisations in each province and serological data. We find the model adequately describes these data, but the addition of interprovincial mobility was not necessary to obtain an accurate description of the 2020-2021 sars-cov-2 pandemic in Belgium. We further demonstrate how our model can be used to help policymakers decide on the optimal timing of the release of social restrictions.We find that adding spatial heterogeneity by geographically stratifying the model results in more uncertain model projections as compared to an equivalent nation-level model, which has both communicative advantages and disadvantages. We finally discuss the impact of imposing local mobility or social contact restrictions to contain an epidemic in a given province and find that lowering social contact is a more effective strategy than lowering mobility.

Estimating the fitness cost and benefit of antimicrobial resistance from pathogen genomic data

Increasing levels of antibiotic resistance in many bacterial pathogen populations are a major threat to public health. Resistance to an antibiotic provides a fitness benefit when the bacteria are exposed to this antibiotic, but resistance also often comes at a cost to the resistant pathogen relative to susceptible counterparts. We lack a good understanding of these benefits and costs of resistance for many bacterial pathogens and antibiotics, but estimating them could lead to better use of antibiotics in a way that reduces or prevents the spread of resistance. Here, we propose a new model for the joint epidemiology of susceptible and resistant variants, which includes explicit parameters for the cost and benefit of resistance. We show how Bayesian inference can be performed under this model using phylogenetic data from susceptible and resistant lineages and that by combining data from both we are able to disentangle and estimate the resistance cost and benefit parameters separately. We applied our inferential methodology to several simulated datasets to demonstrate good scalability and accuracy. We analysed a dataset of Neisseria gonorrhoeae genomes collected between 2000 and 2013 in the USA. We found that two unrelated lineages resistant to fluoroquinolones shared similar epidemic dynamics and resistance parameters. Fluoroquinolones were abandoned for the treatment of gonorrhoea due to increasing levels of resistance, but our results suggest that they could be used to treat a minority of around 10% of cases without causing resistance to grow again.

A tractable physical model for the yeast polarity predicts epistasis and fitness

Accurate phenotype prediction based on genetic information has numerous societal applications, such as crop design or cellular factories. Epistasis, when biological components interact, complicates modelling phenotypes from genotypes. Here we show an approach to mitigate this complication for polarity establishment in budding yeast, where mechanistic information is abundant. We coarse-grain molecular interactions into a so-called mesotype, which we combine with gene expression noise into a physical cell cycle model. First, we show with computer simulations that the mesotype allows validation of the most current biochemical polarity models by quantitatively matching doubling times. Second, the mesotype elucidates epistasis emergence as exemplified by evaluating the predicted mutational effect of key polarity protein Bem1p when combined with known interactors or under different growth conditions. This example also illustrates how unlikely evolutionary trajectories can become more accessible. The tractability of our biophysically justifiable approach inspires a road-map towards bottom-up modelling complementary to statistical inferences.

This article is part of the theme issue ‘Interdisciplinary approaches to predicting evolutionary biology’.

Tumour growth control: analysis of alternative approaches

In this work we address the problem of tumour growth control by properly exploiting a low-dimensional model that grounds on the Chemical Reaction Network (CRN) formalism. Originally conceived to work both in deterministic and stochastic frameworks, it is shown that, except for the case of very low number of tumour cells, the deterministic approach is appropriate to characterize the system behaviour, especially for control planning purposes. Two alternative control approaches are here investigated. One trivially assumes a constant infusion of external drug administration, the other is designed according to a state-feedback control scheme, with complete or partial knowledge of the state. Pros and cons of both control laws are investigated, showing that the tumour size at the beginning of the therapy plays a role of paramount importance for fixed infusion therapies, whilst only state-feedback laws can eradicate arbitrarily large tumours.

The effect of bottleneck size on evolution in nested Darwinian populations

Previous work has shown how a minimal ecological structure consisting of patchily distributed resources and recurrent dispersal between patches can scaffold Darwinian properties onto collections of cells. When the timescale of dispersal is long compared with the time to consume resources, patch fitness increases but comes at a cost to cell growth rates. This creates conditions that initiate evolutionary transitions in individuality. A key feature of the scaffold is a bottleneck created during dispersal, causing patches to be founded by single cells. The bottleneck decreases competition within patches and, hence, creates a strong hereditary link at the level of patches. Here, we construct a fully stochastic model to investigate the effect of bottleneck size on the evolutionary dynamics of both cells and collectives. We show that larger bottlenecks simply slow the dynamics, but, at some point, which depends on the parameters of the within-patch model, the direction of evolution towards the equilibrium reverses. Introduction of random fluctuations in bottleneck sizes with some positive probability of smaller sizes counteracts this, even when the probability of smaller bottlenecks is minimal.

Inferring density-dependent population dynamics mechanisms through rate disambiguation for logistic birth-death processes

Density dependence is important in the ecology and evolution of microbial and cancer cells. Typically, we can only measure net growth rates, but the underlying density-dependent mechanisms that give rise to the observed dynamics can manifest in birth processes, death processes, or both. Therefore, we utilize the mean and variance of cell number fluctuations to separately identify birth and death rates from time series that follow stochastic birth-death processes with logistic growth. Our nonparametric method provides a novel perspective on stochastic parameter identifiability, which we validate by analyzing the accuracy in terms of the discretization bin size. We apply our method to the scenario where a homogeneous cell population goes through three stages: (1) grows naturally to its carrying capacity, (2) is treated with a drug that reduces its carrying capacity, and (3) overcomes the drug effect to restore its original carrying capacity. In each stage, we disambiguate whether the dynamics occur through the birth process, death process, or some combination of the two, which contributes to understanding drug resistance mechanisms. In the case of limited sample sizes, we provide an alternative method based on maximum likelihood and solve a constrained nonlinear optimization problem to identify the most likely density dependence parameter for a given cell number time series. Our methods can be applied to other biological systems at different scales to disambiguate density-dependent mechanisms underlying the same net growth rate.