Journal of Mathematical Biology

This paper concerns the inverse problem of characterising the state of a bioreactor from observations. In laboratory settings, the bioreactor is represented by a device called a chemostat. We consider a differential description of the evolution of the state of the chemostat under environmental fluctuations. First, we model the state evolution as a stochastic process driven by Brownian motion. Under this model, our best knowledge about the state of the chemostat is described by its probability distribution in time, given the distribution of the initial state. The corresponding probability density function solves a deterministic partial differential equation (PDE), the Kolmogorov forward equation. While this provides a probabilistic description, incorporating an observation process allows for a more refined characterisation of the state. More formally, we are interested in obtaining the distribution of the state conditional on an observation process as the solution to a filtering problem, with the corresponding conditional probability density function solving a non-linear stochastic PDE, the Kushner-Stratonovich equation. This paper focuses on the pathwise formulation of this filtering problem in which inferences about the state are obtained conditional on a fixed stream of observations. We establish the existence and uniqueness of solutions to the governing differential equations, ensuring well-posedness before presenting numerical approximations. We approximate the pathwise solution to the filtering problem by combining the finite difference and splitting methods for solving PDEs, and then compare the approximated solution with results from a linearisation method and a classical sequential Monte Carlo method.

Since its introduction by Briat, Gupta and Khammash, the antithetic feedback controller design has attracted considerable attention in both theoretical and experimental systems biology. The case in which the plant is a two-dimensional linear system (making the closed-loop system a nonlinear four-dimensional system) has been analyzed in much detail. This system has a unique equilibrium but, depending on parameters, it may exhibit periodic orbits. An interesting open question is whether other dynamical behaviors, such as chaotic attractors, might be possible for some parameter choices. This note shows that, for any parameter choices, every bounded trajectory satisfies a Poincare-Bendixson property. The analysis is based on the recently introduced notion of k-cooperative dynamical systems. It is shown that the model is a strongly 2-cooperative system, implying that the dynamics in the omega-limit set of any precompact solution is conjugate to the dynamics in a compact invariant subset of a two-dimensional Lipschitz dynamical system, thus precluding chaotic and other strange attractors.

A fresh approach to the dynamics of gene assemblies is presented. Central to the exposition are the concepts of: high value genes; correlated activity; and the orderly unfolding of gene dynamics; and especially dynamic mode decomposition, DMD, a remarkable new tool for dissecting dynamics. This program is carried out, in detail, for the Orlando et al yeast database (Orlando et al. 2008). It is shown that the yeast cell division cycle, CDC, requires no more than a six dimensional space, formed by three complex temporal modal pairs, each associated with characteristic aspects of the cell cycle: (1) A mother cell cohort that follows a fast clock; (2) A daughter cell cohort that follows a slower clock; (3) inherent gene expression, unrelated to the CDC. A derived set of sixty high-value genes serves as a model for the correlated unfolding of gene activity. Confirmation of our results comes from an independent database, and other considerations. The present analysis, leads naturally, to a Fourier description, for the sparsely sampled data. From this, resolved peak times of gene expression are obtained. This in turn leads to prediction of precise times of expression in the unfolding of the CDC genes. The activation of each gene appears as uncoupled dynamics from the mother and daughter cohorts, of different durations. These deliberations lead to detailed estimates of the fraction of mother and daughter cells, specific estimates of their maturation periods, and specific estimates of the number of genes in these cells. An algorithmic framework for yeast modeling is proposed, and based on the new analyses, a range of theoretical ideas and new experiments are suggested. A Supplement contains additional material and other perspectives.

We introduce a mechanistic, nonlocal tumour-growth model designed specifically to capture explosive dynamics that are not adequately explained by standard logistic reaction-diffusion descriptions. The motivation is empirical: the universal scaling law reported in [1] provides compelling cross-sectional evidence of superlinear tumour activity versus tumour burden, but as a phenomenological relationship it does not by itself supply a dynamical mechanism, nor does it rigorously describe how explosive growth emerges, how fast it develops, or how spatial interactions and tissue boundaries influence it. Our model addresses this gap by incorporating nonlocal proliferative feedback--cells respond to a spatially aggregated neighbourhood signal--and a singular, Kawarada-type acceleration that produces "quenching": tumour density stays bounded while the proliferative drive becomes unbounded as the aggregated signal approaches a critical threshold. This offers a concrete mechanistic route to explosive escalation consistent with physical boundedness. We analyse the model under no-flux (Neumann) boundary conditions, appropriate for reflecting tissue interfaces. In the spatially homogeneous setting we prove finite-time onset of the explosive regime and obtain explicit rates for how rapidly it is approached. For spatially heterogeneous perturbations we derive a transparent spectral stability theory showing how the interaction kernel selects spatial scales and how the singular acceleration tightens stability margins as the explosive threshold is approached. These results provide interpretable links between nonlocal interaction structure, boundary effects, and the emergence of rapid growth. Finally, to connect mechanism to data in the spirit of [1], we embed the model in a Bayesian inference framework that treats the interaction kernel and the acceleration strength as unknown and learned from tumour-growth observations. This enables uncertainty-aware estimation of explosive onset times, escalation rates, and stability margins, while positioning the scaling law of [1] as an observable signature that our mechanistic model can explain and quantify rather than merely fit.

The exact mechanisms behind many morphogenic processes are still a mystery. Mechanical cues, such as curvature, play an important role when tissue or cell shape is formed. In this work, we derive and analyze a mechanochemical model. This particular spatially one-dimensional model describes the deformation of a tissue- or cell surface over time, which is driven by a morphogen that locally induces curvature. The model consists of two PDEs with periodic boundary conditions; one reaction-diffusion equation for the morphogen and one PDE that describes the dynamics of the curve, derived by taking the L2-gradient flow of the Helfrich energy. We analyze the possible steady states of this model using geometric singular perturbation theory. It turns out that the strength of interaction between the morphogen and the curvature plays a key role in the type of possible steady state solutions. In the case of weak interaction, the geometry of the slow manifolds allows only for (in space) slowly changing periodic orbits that lay completely on one slow manifold. In the case of strong interaction, there exist multiple front solutions: periodic orbits that jump between different slow manifolds. The singular skeletons of the steady state solutions do not meet the required consistency conditions for the curvature, a priori indicating that the solutions might not be observable. The observability and stability are investigated further using numerical simulation.

It is well known that the presence of an incoherent feedforward loop (IFFL) in a network may give rise to a steady state non-monotonic dose response. This note shows that the converse implication does not hold. It gives an example of a three-dimensional system that has no IFFLs, yet its dose response is bell-shaped. It also studies under what conditions the result is true for two-dimensional systems, in the process recovering, in far more generality, a result given in the T-cell activation literature.

In tissue engineering, it is important to conceive and construct artificial bio-mimetic scaffolds able to foster cell migration as this is a fundamental process in wound healing and tissue regeneration. In order to do that, cubically symmetric and triply periodic porous structures have been identified as promising candidates for instance for the reconstruction of artificial cartilages and bones, also due to their tunable mechanical characteristics and highly inter-connected porous architectures that mimic the trabecular bone hyperboloidal topography. We propose here a mathematical approach that might be helpful to identify what are the best geometrical characteristics of such scaffolds, in order to promote cell migration into the porous structures and speed-up their re-population. The method is based on the observation that cell nucleus deformations should be avoided, yet assuring a good possibility for the cell to reach the wall of the porous structure. Mathematically speaking, this leads to the problem of identifying the size of the largest sphere that can pass, without being stuck, through the pores of the bio-mimetic scaffold.

The analysis of large gene and metabolic networks is often hindered by unknown biochemical parameters and the nonlinear nature of classical S-system models. To address this, we introduce a framework based on combinatorial toric geometry computed with tools such as Normaliz, SageMath, it is worth mentioning this technique in not restrictive to integer vectors, there exists a natural extension to real geometries. Unlike traditional approaches, which rely on parameter dependent fixed points, our method constructs a Topological Environment derived from the dual space of kinetic orders, leading to what we call orthogonal enzyme kinetics. Within this topological setting, fixed points are computed on the algebraic torus, enabling the transformation of nonlinear dynamics into linear forms. Importantly, these fixed points are independent of kinetic parameters and depend only on network topology and interaction signs. Applying this methodology to gene circuits involved in circadian rhythms, we reproduce previously reported oscillatory physiologies.

This note describes the outcome of an epidemic in a heterogeneous population with a very simple structure. The population is split into two sections in which the epidemic runs a different course. The reproduction numbers of the two epidemics are unobserved, only the overall reproduction number is known. For such a population the outcome of the epidemic can be as well far worse as far better than expected on the basis of the overall reproduction number. By considering a very simple model of this population some calculations are feasible under general assumptions on the epidemic itself. These calculations show in which direction models based on the overall reproduction number can misrepresent the real-world situation.

We discuss some critical events of the origins of life using a mathematical model and simulation studies. We find that for a replicating population of RNA molecules participating in template-directed polymerization, the hitting and establishment of a high-fidelity replicator depends critically on the polymerase fitness and sequence specificity landscapes and on genome dimension. Probability of hitting is dominated by polymerase landscape curvature, whereas hitting time is dominated by genome dimension. Surface chemistries, compartmentalization, and decay increase hitting times. These results suggest replication to be the first privileged function marking the start of Darwinian evolution, possibly in conjunction with clay minerals or preceded by metabolism, whose dynamics evolved mostly during the final period of the search.

Cancer progression is driven by the interplay between cancer cell phenotypic variability regulated by EMT/MET, and cell-cell and cell-matrix interactions. Starting with an individual-based model, in which every cell is characterised by its position, velocity and a continuously varying epithelial-mesenchymal phenotype, we derive, via a kinetic description and a mean-field limit, two alternative macroscopic formulations: an Euler-like system that couples mass and momentum to the phenotypic variable, and a single advection-aggregation-diffusion equation (AADE) for the cancer cell density. Both macroscopic models retain the non-local adhesion-repulsion forces, haptotactic response to an evolving ECM, and phenotype-dependent transition dynamics driven by TGF-{beta}. Numerical experiments in two spatial dimensions indicate that the macroscopic equations reproduce key scenarios obtained at the individual scale. In particular, a microscopic-macroscopic comparison shows that the AADE density reproduces acurately both the spatial localisation and the phenotypic decomposition of the individual cell population. We also demonstrate that varying only the steepness of the TGF-{beta} switch function, changes the EMT response from an almost binary epithelial-mesenchymal (EM) separation to a partial EM phenotypes. This study provides a systematic bridge from stochastic, heterogeneous cell dynamics to continuum descriptions, for investigating phenotype driven tumour invasion and supporting the choice of macroscopic models in large-scale simulations and analytical studies.

The global potential of a chemical reaction network has many applications and is closely related to the stochastic detailed balance. However, many fundamental questions concerning stochastic detailed balance remain unresolved, such as whether it depends on the system volume and how to construct new systems that satisfy it. In this paper, we show that stochastic detailed balance may depend on the system volume. We therefore introduce four types of stochastic detailed balance according to their dependence on volume and rate constants, and systematically investigate the relationships among them. Our results distinguish detailed balance arising from particular choices of volume and parameters from that enforced by network structure, and identify conditions under which detailed balance at one volume extends to all volumes. We further obtain a class of networks satisfying stochastic detailed balance for every volume and every positive choice of rate constants, and construct new systems whose global potentials exhibit double-well structures.

The tetrapod limb skeleton is initiated in unpatterned limb bud mesenchyme by the formation of precartilage condensations. Here, based on time-lapse videographic analysis of a forming condensation in a high-density culture of chicken limb bud mesenchyme, we observe a phase transition to a more fluidized state for cells within spatial compacted foci (protocondensations that will progress to condensations), as reflected in their spatial confinement, cell-substratum interaction and speed of motion. Previous work showed that galectin-8 and galectin-1A, two proteins of the galactoside-binding galectin family, are the earliest determinants of this process in the chicken limb bud, and that their interactions in forming skeletogenic patterns of condensations can be interpreted mathematically through a reaction-diffusion-adhesion framework. Based on this framework, we use an ordinary differential equation-based approach to analyze the core switching modality of the galectin reaction network and characterize the states of the network independent of the diffusive and adhesive arms of the patterning mechanism. We identify two steady states where the concentrations of both galectins are respectively, negligible, and very high. An explicit Lyapunov function shows that there are no periodic solutions. For sigmoidal galectin production terms, the model exhibits a bistable switch that arises from a monostable state via saddle-node bifurcation. Our model therefore predicts that the galectin network exists in low and high expression states separated in space or time without any intermediate states. This provides a causal basis for the observed outside vs. inside transition observed in the in vitro video data. We performed a quantitative analysis of the distribution of galectin-1A in cultures of condensing chick limb mesenchymal cells and found that the interior of the protocondensations had concentrations of this protein (compared to the immediate exterior) over and above that expected from its higher cell density, consistent with the models predictions. The galectin-based patterning network is thus suggested, on theoretical grounds, to incorporate a core switch independent of any spatial or temporal dynamics, that drives the chondrogenic cell state transition in limb skeletogenesis.

This paper discusses the mean duration of a closed epidemic modeled by a discrete-time Markov chain. We develop a methodology for the efficient computation of the quantity of interest. The Markov chain model in consideration is bivariate, and is formally handled. We derive explicit terms for the probability to transition from one state to another, and prove that the chain is absorbing. The computation of the mean duration is translated to the computation of the expected hitting times to the set of absorbing states. We use the theory of absorbing Markov chains to derive a matrix formulation that gives way to an efficient algorithm to solve for the expected hitting times. This approach is instantiated in the form of a concrete algorithm, which is further optimized by using dynamic programming. Finally, we have implemented the method and tested it against the use of simulations to estimate mean durations.

A survey of the literature reveals notable discrepancies among the purported exact results for the spectra of stochastic gene expression models. For self-repressing gene circuits, previous studies ([Phys. Rev. Lett. 99, 108103 (2007)], [Phys. Rev. E 83,062902 (2011)], [J. Chem. Phys. 160, 074105 (2024)], and [bioRxiv 2025.02.05.635946 (2025)]) have provided different exact solutions for the eigenvalues of the generator matrix. In this work, we propose a unified Hilbert space framework for the spectral theory of stochastic gene expression. Based on this framework, we analytically derive the spectra for models of constitutive, bursty, and autoregulated gene expression. The eigenvalues and eigenvectors obtained are then used to construct an exact spectral representation of the time-dependent distribution of gene product numbers. The spectral gap between the zero eigenvalue and the first nonzero eigenvalue, which reflects the relaxation rate of the system towards its steady state, is then compared with the prediction of the deterministic model, and we find that deterministic modeling fails to capture the relaxation rate when autoregulation is strong. In particular, our results demonstrate that for infinite-dimensional operators such as in stochastic gene expression models, many conclusions in linear algebra do not apply, and one must rely on the modern theory of functional analysis.

Surprisingly, the discrete-time version of the general 1927 Kermack-McKendrick epidemic model has, to our knowledge, not been formulated in the literature, and we rectify this omission here. The discrete time version is as general and flexible as its continuous-time counterpart, and contains numerous compartmental models as special cases. In contrast to the continuous time version, the discrete time version of the model is very easy to implement computationally, and thus promises to become a powerful tool for exploring control scenarios for specific infectious diseases. To demonstrate the potential, we investigate numerically how the incidence-peak size depends on model ingredients. We find that, with the same reproduction number and initial speed of epidemic spread, compartmental models systematically predict lower peak sizes than models that use a fixed duration for the latent and infectious periods.

Melanoma is a cancer of the melanocyte, known to have an ability to readily switch between different transcriptional cell states that convey different phenotypic properties (e.g. hyper-differentiated, neural crest-like). This ability is believed to underpin intratumour heterogeneity and plastic adaptation, which contributes to resistance to therapy and immune evasion of the tumour. Therefore, understanding the mechanisms underlying acquisition of transcriptional cell states and cell-state switching is crucial for the development of therapies. We model a minimal gene regulatory network comprising three key transcription factors, whose varying gene expression encodes different melanoma cell states, and use deterministic spatiotemporal differential-equation models to study gene-expression dynamics. We exploit an approximation, based on cooperative binding of transcription factors, in which the models are piecewise-linear. We classify stable states of the local model in a biologically relevant manner and, using a naive model of intercellular communication, we explore how a population of cells can take on a shared characteristic through travelling waves of gene expression. We derive a condition determining which characteristic will become dominant, under sufficiently strong cell-cell signalling, which creates a partition of parameter space.

AO_SCPLOWBSTRACTC_SCPLOWCompartmental models have long served as important tools in mathematical epidemiology, with their usefulness highlighted by the recent COVID-19 pandemic. However, most of the classical models fail to account for certain features of this disease and others like it, such as the ability of exposed individuals to recover without becoming infectious, or the possibility that asymptomatic individuals can indeed transmit the disease but at a lesser rate than the symptomatic. Furthermore, the rise of new disease variants and the imperfection of vaccines suggest that concept of endemic equilibrium is perhaps more pertinent than that of herd immunity. Here we propose a new compartmental epidemiological model and study its equilibria, characterizing the stability of both the endemic and disease-free equilibria in terms of the basic reproductive number. Moreover, we introduce a second compartmental model, generalizing our first, which accounts for vaccinated individuals, and begin an analysis of its equilibria.

In a variety of practical applications, there is a need to investigate diffusion or reaction-diffusion processes on complex structures, including brain networks, that can be modeled as weighted undirected and directed graphs. As an instance, the celebrated Fisher-Kolmogorov-Petrovsky-Piskunov (Fisher-KPP) reaction-diffusion equation are becoming increasingly popular for use in graph frameworks by substituting the standard graph Laplacian operator for the continuous one to study the progression of neurodegenerative diseases such as tauopathies including Alzheimers disease (AD). However, due to the porous structure of neuronal fibers, the spreading of toxic species can be governed by an anomalous diffusion process rather than a normal one, and if this is the case, the standard graph Laplacian cannot adequately describe the dynamics of the spreading process. To capture such more complicated dynamics, we propose a diffusion equation with a nonlinear Laplacian operator and a generalization of the Fisher-KPP reaction-diffusion equation on undirected and directed networks using extensions of fractional polynomial (FP) functions. A complete analysis is also provided for the extended FP diffusion equation, including existence, uniqueness, and convergence of solutions, as well as stability of equilibria. Moreover, for the extended FP Fisher-KPP reaction-diffusion equation, we derive a family of positively invariant sets allowing us to establish existence, uniqueness, and boundedness of solutions. Finally, we conclude by investigating nonlinear diffusion on a directed one-dimensional lattice and then modeling tauopathy progression in the mouse brain to gain a deeper understanding of the potential applications of the proposed extended FP equations.

Matrix games under time constraints are natural extensions of matrix games. They consider the fact that, in addition to the payoff, a pairwise interaction has a further consequence for the contestants. Namely, both players have to wait for some time before becoming fit to participate in a subsequent interaction. Every matrix game can be assigned a continuous dynamical system (the replicator equation) which describes how the frequencies of different phenotypes evolve in the population. One of the fundamental theorems of evolutionary matrix games asserts that the state corresponding to an evolutionarily stable strategy is an asymptotically stable rest point of the replicator equation (Taylor and Yonker 1978, Hofbauer et al. 1979, Zeeman 1980). Garay et al. (2018) and Varga et al. (2020) generalized the statement to two-strategy and, in some particular cases, three- or more strategy matrix games under time constraints. However, the question of whether the implication holds in general remained open. Here examples are provided demonstrating that the answer is no. Moreover, we point out through the rock-scissor-paper game that arbitrary small differences between waiting times can destabilize the rest point corresponding to an ESS. It is also shown that a stable limit cycle can arise around the unstable rest point in a supercritical Hopf bifurcation. Mathematics Subject Classification91A22, 92D15, 92D25, 91A80, 91A05, 91A10, 91A40, 92D40