A two-strain model of infectious disease spread with asymmetric temporary immunity periods and partial cross-immunity.

We introduce a two-strain model with asymmetric temporary immunity periods and partial cross-immunity. We derive explicit conditions for competitive exclusion and coexistence of the strains depending on the strain-specific basic reproduction numbers, temporary immunity periods, and degree of cross-immunity. The results of our bifurcation analysis suggest that, even when two strains share similar basic reproduction numbers and other epidemiological parameters, a disparity in temporary immunity periods and partial or complete cross-immunity can provide a significant competitive advantage. To analyze the dynamics, we introduce a quasi-steady state reduced model which assumes the original strain remains at its endemic steady state. We completely analyze the resulting reduced planar hybrid switching system using linear stability analysis, planar phase-plane analysis, and the Bendixson-Dulac criterion. We validate both the full and reduced models with COVID-19 incidence data, focusing on the Delta (B.1.617.2), Omicron (B.1.1.529), and Kraken (XBB.1.5) variants. These numerical studies suggest that, while early novel strains of COVID-19 had a tendency toward dramatic takeovers and extinction of ancestral strains, more recent strains have the capacity for co-existence.

A framework for deriving analytic steady states of biochemical reaction networks.

The long-term behaviors of biochemical systems are often described by their steady states. Deriving these states directly for complex networks arising from real-world applications, however, is often challenging. Recent work has consequently focused on network-based approaches. Specifically, biochemical reaction networks are transformed into weakly reversible and deficiency zero generalized networks, which allows the derivation of their analytic steady states. Identifying this transformation, however, can be challenging for large and complex networks. In this paper, we address this difficulty by breaking the complex network into smaller independent subnetworks and then transforming the subnetworks to derive the analytic steady states of each subnetwork. We show that stitching these solutions together leads to the analytic steady states of the original network. To facilitate this process, we develop a user-friendly and publicly available package, COMPILES (COMPutIng anaLytic stEady States). With COMPILES, we can easily test the presence of bistability of a CRISPRi toggle switch model, which was previously investigated via tremendous number of numerical simulations and within a limited range of parameters. Furthermore, COMPILES can be used to identify absolute concentration robustness (ACR), the property of a system that maintains the concentration of particular species at a steady state regardless of any initial concentrations. Specifically, our approach completely identifies all the species with and without ACR in a complex insulin model. Our method provides an effective approach to analyzing and understanding complex biochemical systems.

A data-validated temporary immunity model of COVID-19 spread in Michigan.

We introduce a distributed-delay differential equation disease spread model for COVID-19 spread. The model explicitly incorporates the population's time-dependent vaccine uptake and incorporates a gamma-distributed temporary immunity period for both vaccination and previous infection. We validate the model on COVID-19 cases and deaths data from the state of Michigan and use the calibrated model to forecast the spread and impact of the disease under a variety of realistic booster vaccine strategies. The model suggests that the mean immunity duration for individuals after vaccination is 350 days and after a prior infection is 242 days. Simulations suggest that both high population-wide adherence to vaccination mandates and a more-than-annually frequency of booster doses will be required to contain outbreaks in the future.

A dynamical framework for modeling fear of infection and frustration with social distancing in COVID-19 spread.

We introduce a novel modeling framework for incorporating fear of infection and frustration with social distancing into disease dynamics. We show that the resulting SEIR behavior-perception model has three principal modes of qualitative behavior-no outbreak, controlled outbreak, and uncontrolled outbreak. We also demonstrate that the model can produce transient and sustained waves of infection consistent with secondary outbreaks. We fit the model to cumulative COVID-19 case and mortality data from several regions. Our analysis suggests that regions which experience a significant decline after the first wave of infection, such as Canada and Israel, are more likely to contain secondary waves of infection, whereas regions which only achieve moderate success in mitigating the disease's spread initially, such as the United States, are likely to experience substantial secondary waves or uncontrolled outbreaks.

Realizations of kinetic differential equations.

The induced kinetic differential equations of a reaction network endowed with mass action type kinetics is a system of polynomial differential equations. The problem studied here is: Given a system of polynomial differential equations, is it possible to find a network which induces these equations; in other words: is it possible to find a kinetic realization of this system of differential equations? If yes, can we find a network with some chemically relevant properties (implying also important dynamic consequences), such as reversibility, weak reversibility, zero deficiency, detailed balancing, complex balancing, mass conservation, etc.? The constructive answers presented to a series of questions of the above type are useful when fitting differential equations to datasets, or when trying to find out the dynamic behavior of the solutions of differential equations. It turns out that some of these results can be applied when trying to solve seemingly unrelated mathematical problems, like the existence of positive solutions to algebraic equations.

Computing Weakly Reversible Deficiency Zero Network Translations Using Elementary Flux Modes.

We present a computational method for performing structural translation, which has been studied recently in the context of analyzing the steady states and dynamical behavior of mass-action systems derived from biochemical reaction networks. Our procedure involves solving a binary linear programming problem where the decision variables correspond to interactions between the reactions of the original network. We call the resulting network a reaction-to-reaction graph and formalize how such a construction relates to the original reaction network and the structural translation. We demonstrate the efficacy and efficiency of the algorithm by running it on 508 networks from the European Bioinformatics Institutes' BioModels database. We also summarize how this work can be incorporated into recently proposed algorithms for establishing mono- and multistationarity in biochemical reaction systems.

A Deficiency-Based Approach to Parametrizing Positive Equilibria of Biochemical Reaction Systems.

We present conditions which guarantee a parametrization of the set of positive equilibria of a generalized mass-action system. Our main results state that (1) if the underlying generalized chemical reaction network has an effective deficiency of zero, then the set of positive equilibria coincides with the parametrized set of complex-balanced equilibria and (2) if the network is weakly reversible and has a kinetic deficiency of zero, then the equilibrium set is nonempty and has a positive, typically rational, parametrization. Via the method of network translation, we apply our results to classical mass-action systems studied in the biochemical literature, including the EnvZ-OmpR and shuttled WNT signaling pathways. A parametrization of the set of positive equilibria of a (generalized) mass-action system is often a prerequisite for the study of multistationarity and allows an easy check for the occurrence of absolute concentration robustness, as we demonstrate for the EnvZ-OmpR pathway.

Network Translation and Steady-State Properties of Chemical Reaction Systems.

Network translation has recently been used to establish steady-state properties of mass action systems by corresponding the given system to a generalized one which is either dynamically or steady-state equivalent. In this work, we further use network translation to identify network structures which give rise to the well-studied property of absolute concentration robustness in the corresponding mass action systems. In addition to establishing the capacity for absolute concentration robustness, we show that network translation can often provide a method for deriving the steady-state value of the robust species. We furthermore present a MILP algorithm for the identification of translated chemical reaction networks that improves on previous approaches, allowing for easier application of the theory.

Conditions for extinction events in chemical reaction networks with discrete state spaces.

We study chemical reaction networks with discrete state spaces and present sufficient conditions on the structure of the network that guarantee the system exhibits an extinction event. The conditions we derive involve creating a modified chemical reaction network called a domination-expanded reaction network and then checking properties of this network. Unlike previous results, our analysis allows algorithmic implementation via systems of equalities and inequalities and suggests sequences of reactions which may lead to extinction events. We apply the results to several networks including an EnvZ-OmpR signaling pathway in Escherichia coli.

A computational approach to extinction events in chemical reaction networks with discrete state spaces.

Recent work of Johnston et al. has produced sufficient conditions on the structure of a chemical reaction network which guarantee that the corresponding discrete state space system exhibits an extinction event. The conditions consist of a series of systems of equalities and inequalities on the edges of a modified reaction network called a domination-expanded reaction network. In this paper, we present a computational implementation of these conditions written in Python and apply the program on examples drawn from the biochemical literature. We also run the program on 458 models from the European Bioinformatics Institute's BioModels Database and report our results.

A computational approach to persistence, permanence, and endotacticity of biochemical reaction systems.

We introduce a mixed-integer linear programming (MILP) framework capable of determining whether a chemical reaction network possesses the property of being endotactic or strongly endotactic. The network property of being strongly endotactic is known to lead to persistence and permanence of chemical species under genetic kinetic assumptions, while the same result is conjectured but as yet unproved for general endotactic networks. The algorithms we present are the first capable of verifying endotacticity of chemical reaction networks for systems with greater than two constituent species. We implement the algorithms in the open-source online package CoNtRol and apply them to a large sample of networks from the European Bioinformatics Institute's BioModels Database. We use strong endotacticity to establish for the first time the permanence of a well-studied circadian clock mechanism.

A Computational Approach to Steady State Correspondence of Regular and Generalized Mass Action Systems.

It has been recently observed that the dynamical properties of mass action systems arising from many models of biochemical reaction networks can be characterized by considering the corresponding properties of a related generalized mass action system. The correspondence process known as network translation in particular has been shown to be useful in characterizing a system's steady states. In this paper, we further develop the theory of network translation with particular focus on a subclass of translations known as improper translations. For these translations, we derive conditions on the network topology of the translated network which are sufficient to guarantee the original and translated systems share the same steady states. We then present a mixed-integer linear programming algorithm capable of determining whether a mass action system can be corresponded to a generalized system through the process of network translation.

Translated chemical reaction networks.

Many biochemical and industrial applications involve complicated networks of simultaneously occurring chemical reactions. Under the assumption of mass action kinetics, the dynamics of these chemical reaction networks are governed by systems of polynomial ordinary differential equations. The steady states of these mass action systems have been analyzed via a variety of techniques, including stoichiometric network analysis, deficiency theory, and algebraic techniques (e.g., Gröbner bases). In this paper, we present a novel method for characterizing the steady states of mass action systems. Our method explicitly links a network's capacity to permit a particular class of steady states, called toric steady states, to topological properties of a generalized network called a translated chemical reaction network. These networks share their reaction vectors with their source network but are permitted to have different complex stoichiometries and different network topologies. We apply the results to examples drawn from the biochemical literature.

Stochastic analysis of biochemical reaction networks with absolute concentration robustness.

It has recently been shown that structural conditions on the reaction network, rather than a 'fine-tuning' of system parameters, often suffice to impart 'absolute concentration robustness' (ACR) on a wide class of biologically relevant, deterministically modelled mass-action systems. We show here that fundamentally different conclusions about the long-term behaviour of such systems are reached if the systems are instead modelled with stochastic dynamics and a discrete state space. Specifically, we characterize a large class of models that exhibit convergence to a positive robust equilibrium in the deterministic setting, whereas trajectories of the corresponding stochastic models are necessarily absorbed by a set of states that reside on the boundary of the state space, i.e. the system undergoes an extinction event. If the time to extinction is large relative to the relevant timescales of the system, the process will appear to settle down to a stationary distribution long before the inevitable extinction will occur. This quasi-stationary distribution is considered for two systems taken from the literature, and results consistent with ACR are recovered by showing that the quasi-stationary distribution of the robust species approaches a Poisson distribution.

Computing weakly reversible linearly conjugate chemical reaction networks with minimal deficiency.

Mass-action kinetics is frequently used in systems biology to model the behavior of interacting chemical species. Many important dynamical properties are known to hold for such systems if their underlying networks are weakly reversible and have a low deficiency. In particular, the Deficiency Zero and Deficiency One Theorems guarantee strong regularity with regards to the number and stability of positive equilibrium states. It is also known that chemical reaction networks with distinct reaction structure can admit mass-action systems with the same qualitative dynamics. The theory of linear conjugacy encapsulates the cases where this relationship is captured by a linear transformation. In this paper, we propose a mixed-integer linear programming algorithm capable of determining the minimal deficiency weakly reversible reaction network which admits a mass-action system which is linearly conjugate to a given reaction network.

On the proportion of cancer stem cells in a tumour.

It is now generally accepted that cancers contain a sub-population, the cancer stem cells (CSCs), which initiate and drive a tumour's growth. At least until recently it has been widely assumed that only a small proportion of the cells in a tumour are CSCs. Here we use a mathematical model, supported by experimental evidence, to show that such an assumption is unwarranted. We show that CSCs may comprise any possible proportion of the tumour, and that the higher the proportion the more aggressive the tumour is likely to be.

Examples of mathematical modeling: tales from the crypt.

Mathematical modeling is being increasingly recognized within the biomedical sciences as an important tool that can aid the understanding of biological systems. The heavily regulated cell renewal cycle in the colonic crypt provides a good example of how modeling can be used to find out key features of the system kinetics, and help to explain both the breakdown of homeostasis and the initiation of tumorigenesis. We use the cell population model by Johnston et al. to illustrate the power of mathematical modeling by considering two key questions about the cell population dynamics in the colonic crypt. We ask: how can a model describe both homeostasis and unregulated growth in tumorigenesis; and to which parameters in the system is the model most sensitive? In order to address these questions, we discuss what type of modeling approach is most appropriate in the crypt. We use the model to argue why tumorigenesis is observed to occur in stages with long lag phases between periods of rapid growth, and we identify the key parameters.

Mathematical modeling of cell population dynamics in the colonic crypt and in colorectal cancer.

Colorectal cancer is initiated in colonic crypts. A succession of genetic mutations or epigenetic changes can lead to homeostasis in the crypt being overcome, and subsequent unbounded growth. We consider the dynamics of a single colorectal crypt by using a compartmental approach [Tomlinson IPM, Bodmer WF (1995) Proc Natl Acad Sci USA 92:], which accounts for populations of stem cells, differentiated cells, and transit cells. That original model made the simplifying assumptions that each cell population divides synchronously, but we relax these assumptions by adopting an age-structured approach that models asynchronous cell division, and by using a continuum model. We discuss two mechanisms that could regulate the growth of cell numbers and maintain the equilibrium that is normally observed in the crypt. The first will always maintain an equilibrium for all parameter values, whereas the second can allow unbounded proliferation if the net per capita growth rates are large enough. Results show that an increase in cell renewal, which is equivalent to a failure of programmed cell death or of differentiation, can lead to the growth of cancers. The second model can be used to explain the long lag phases in tumor growth, during which new, higher equilibria are reached, before unlimited growth in cell numbers ensues.

Dynamics of vaccination strategies via projected dynamical systems.

Previous game theoretical analyses of vaccinating behaviour have underscored the strategic interaction between individuals attempting to maximise their health states, in situations where an individual's health state depends upon the vaccination decisions of others due to the presence of herd immunity. Here, we extend such analyses by applying the theories of variational inequalities (VI) and projected dynamical systems (PDS) to vaccination games. A PDS provides a dynamics that gives the conditions for existence, uniqueness and stability properties of Nash equilibria. In this paper, it is used to analyse the dynamics of vaccinating behaviour in a population consisting of distinct social groups, where each group has different perceptions of vaccine and disease risks. In particular, we study populations with two groups, where the size of one group is strictly larger than the size of the other group (a majority/minority population). We find that a population with a vaccine-inclined majority group and a vaccine-averse minority group exhibits higher average vaccine coverage than the corresponding homogeneous population, when the vaccine is perceived as being risky relative to the disease. Our model also reproduces a feature of real populations: In certain parameter regimes, it is possible to have a majority group adopting high vaccination rates and simultaneously a vaccine-averse minority group adopting low vaccination rates. Moreover, we find that minority groups will tend to exhibit more extreme changes in vaccinating behaviour for a given change in risk perception, in comparison to majority groups. These results emphasise the important role played by social heterogeneity in vaccination behaviour, while also highlighting the valuable role that can be played by PDS and VI in mathematical epidemiology.