Mathematical model of a cytokine storm.

Cytokine storm is a life-threatening inflammatory response that is characterized by hyperactivation of the immune system, and which can be caused by various therapies, autoimmune conditions, or pathogens, such as respiratory syndrome coronavirus 2 (SARS-CoV-2), which causes coronavirus disease COVID-19. While initial causes of cytokine storms can vary, late-stage clinical manifestations of cytokine storm converge and often overlap, and therefore a better understanding of how normal immune response turns pathological is warranted. Here we propose a theoretical framework, where cytokine storm phenomenology is captured using a conceptual mathematical model, where cytokines can both activate and regulate the immune system. We simulate normal immune response to infection, and through variation of system parameters identify conditions where, within the frameworks of this model, cytokine storm can arise. We demonstrate that cytokine storm is a transitional regime, and identify three main factors that must converge to result in storm-like dynamics, two of which represent individual-specific characteristics, thereby providing a possible explanation for why some people develop CRS, while others may not. We also discuss possible ecological insights into cytokine-immune interactions and provide mathematical analysis for the underlying regimes. We conclude with a discussion of how results of this analysis can be used in future research.

Stable coevolutionary regimes for genetic parasites and their hosts: you must differ to coevolve.

BACKGROUND: Genetic parasites are ubiquitous satellites of cellular life forms most of which host a variety of mobile genetic elements including transposons, plasmids and viruses. Theoretical considerations and computer simulations suggest that emergence of genetic parasites is intrinsic to evolving replicator systems. RESULTS: Using methods of bifurcation analysis, we investigated the stability of simple models of replicator-parasite coevolution in a well-mixed environment. We first analyze what appears to be the simplest imaginable system of this type, one in which the parasite evolves during the replication of the host genome through a minimal mutation that renders the genome of the emerging parasite incapable of producing the replicase but able to recognize and recruit it for its own replication. This model has only trivial or "semi-trivial", parasite-free equilibria: an inefficient parasite is outcompeted by the host and dies off, whereas an efficient one pushes the host out of existence, leading to the collapse of the entire system. We show that stable host-parasite coevolution (a non-trivial equilibrium) is possible in a modified model where the parasite is qualitatively distinct from the host replicator in that the replication of the parasite depends solely on the availability of the host but not on the carrying capacity of the environment. CONCLUSIONS: We analytically determine the conditions for stable coevolution of genetic parasites and their hosts coevolution in simple mathematical models. It is shown that the evolutionary dynamics of a parasite that initially evolves from the host through the loss of the ability to replicate autonomously must substantially differ from that of the host, for a stable host-parasite coevolution regime to be established.

Cancer immunoediting: A process driven by metabolic competition as a predator-prey-shared resource type model.

It is a well-established fact that tumors up-regulate glucose consumption to meet increasing demands for rapidly available energy by upregulating a purely glycolytic mode of glucose metabolism. What is often neglected is that activated cytotoxic cells of the immune system, integral players in the carcinogenesis process, also come to rely on glycolysis as a primary mode of glucose metabolism. Moreover, while cancer cells can revert back to aerobic metabolism, rapidly proliferating cytotoxic lymphocytes are incapable of performing their function when adequate resources are lacking. Consequently, it is likely that in the tumor microenvironment there may exist competition for shared resources between cancer cells and the cells of the immune system, which may underlie much of tumor-immune dynamics. Proposed here is a model of tumor-immune-glucose interactions, formulated as a predator-prey-common resource type system. The outcome of these interactions ranges from tumor elimination, to tumor dormancy, to unrestrained tumor growth. It is also predicted that the process of tumor escape can be preceded by periods of oscillatory tumor growth. A detailed bifurcation analysis of three subsystems of the model suggest that oscillatory regimes are a result of competition for shared resource (glucose) between the predator (immune cells) and the prey (cancer cells). Existence of competition for nutrients between cancer and immune cells may provide additional mechanistic insight as to why the efficacy of many immunotherapies may be limited.

Pseudo-chaotic oscillations in CRISPR-virus coevolution predicted by bifurcation analysis.

BACKGROUND: The CRISPR-Cas systems of adaptive antivirus immunity are present in most archaea and many bacteria, and provide resistance to specific viruses or plasmids by inserting fragments of foreign DNA into the host genome and then utilizing transcripts of these spacers to inactivate the cognate foreign genome. The recent development of powerful genome engineering tools on the basis of CRISPR-Cas has sharply increased the interest in the diversity and evolution of these systems. Comparative genomic data indicate that during evolution of prokaryotes CRISPR-Cas loci are lost and acquired via horizontal gene transfer at high rates. Mathematical modeling and initial experimental studies of CRISPR-carrying microbes and viruses reveal complex coevolutionary dynamics. RESULTS: We performed a bifurcation analysis of models of coevolution of viruses and microbial host that possess CRISPR-Cas hereditary adaptive immunity systems. The analyzed Malthusian and logistic models display complex, and in particular, quasi-chaotic oscillation regimes that have not been previously observed experimentally or in agent-based models of the CRISPR-mediated immunity. The key factors for the appearance of the quasi-chaotic oscillations are the non-linear dependence of the host immunity on the virus load and the partitioning of the hosts into the immune and susceptible populations, so that the system consists of three components. CONCLUSIONS: Bifurcation analysis of CRISPR-host coevolution model predicts complex regimes including quasi-chaotic oscillations. The quasi-chaotic regimes of virus-host coevolution are likely to be biologically relevant given the evolutionary instability of the CRISPR-Cas loci revealed by comparative genomics. The results of this analysis might have implications beyond the CRISPR-Cas systems, i.e. could describe the behavior of any adaptive immunity system with a heritable component, be it genetic or epigenetic. These predictions are experimentally testable. REVIEWERS' REPORTS: This manuscript was reviewed by Sandor Pongor, Sergei Maslov and Marek Kimmel. For the complete reports, go to the Reviewers' Reports section.

Transitional regimes as early warning signals in resource dependent competition models.

In this paper a question of "how much overconsumption a renewable resource can tolerate" is addressed using a mathematical model, where individuals in a parametrically heterogeneous population not only compete for the common resource but can also contribute to its restoration. Through bifurcation analysis a threshold of system resistance to over-consumers (individuals that take more than they restore) was identified, as well as a series of transitional regimes that the population goes through before it exhausts the common resource and thus goes extinct itself, a phenomenon known as "the tragedy of the commons". It was also observed that (1) for some parameter domains a population can survive or go extinct depending on its initial conditions, (2) under the same set of initial conditions, a heterogeneous population survives longer than a homogeneous population and (3) when the natural decay rate of the common resource is high enough, the population can endure the presence of more aggressive over-consumers without going extinct.

On the asymptotic behaviour of the solutions to the replicator equation.

Selection systems and the corresponding replicator equations model the evolution of replicators with a high level of abstraction. In this paper, we apply novel methods of analysis of selection systems to the replicator equations. To be suitable for the suggested algorithm, the interaction matrix of the replicator equation should be transformed; in particular, the standard singular value decomposition allows us to rewrite the replicator equation in a convenient form. The original n-dimensional problem is reduced to the analysis of asymptotic behaviour of the solutions to the so-called escort system, which in some important cases can be of significantly smaller dimension than the original system. The Newton diagram methods are applied to study the asymptotic behaviour of the solutions to the escort system, when interaction matrix has Rank 1 or 2. A general replicator equation with the interaction matrix of Rank 1 is fully analysed; the conditions are provided when the asymptotic state is a polymorphic equilibrium. As an example of the system with the interaction matrix of Rank 2, we consider the problem from Adams & Sornborger (2007, Analysis of a certain class of replicator equations. J. Math. Biol., 54, 357-384), for which we show, for an arbitrary dimension of the system and under some suitable conditions, that generically one globally stable equilibrium exits on the 1-skeleton of the simplex.

Dynamics of population communities with prey migrations and Allee effects: a bifurcation approach.

The population dynamics of predator-prey systems in the presence of patch-specific predators are explored in a setting where the prey population has access to both habitats. The emphasis is in situations where patch-prey abundance drives prey dispersal between patches, with the fragile prey populations, i.e. populations subject to the Allee effect. The resulting 3D and 4D non-linear systems depending on some parameters, which reflect 'measures' of factors under consideration, support rich dynamics and in particular a diverse number of predator-prey life history outcomes. The model's mathematical analysis is carried out via submodels that focus in lower-dimensional settings. The outcomes depend on and, in fact, are quite sensitive to the structure of the system, the range of parameter values and initial conditions. We show that the system can support multistability and a diverse set of predator-prey life-history dynamics that include rather complex dynamical system outcomes. It is argued that, in general, evolution should favour heterogeneous settings including Allee effects, prey refuges and patch-specific predators.

Myeloid cells in tumour-immune interactions.

Despite highly developed specific immune responses, tumour cells often manage to escape recognition by the immune system, continuing to grow uncontrollably. Experimental work suggests that mature myeloid cells may be central to the activation of the specific immune response. Recognition and subsequent control of tumour growth by the cells of the specific immune response depend on the balance between immature (ImC) and mature (MmC) myeloid cells in the body. However, tumour cells produce cytokines that inhibit ImC maturation, altering the balance between ImC and MmC. Hence, the focus of this manuscript is on the study of the potential role of this inhibiting mechanism on tumour growth dynamics. A conceptual predator-prey type model that incorporates the dynamics and interactions of tumour cells, CD8(+) T cells, ImC and MmC is proposed in order to address the role of this mechanism. The prey (tumour) has a defence mechanism (blocking the maturation of ImC) that prevents the predator (immune system) from recognizing it. The model, a four-dimensional nonlinear system of ordinary differential equations, is reduced to a two-dimensional system using time-scale arguments that are tied to the maturation rate of ImC. Analysis shows that the model is capable of supporting biologically reasonable patterns of behaviour depending on the initial conditions. A range of parameters, where healing without external influences can occur, is identified both qualitatively and quantitatively.

"Traveling wave" solutions of FitzHugh model with cross-diffusion.

The FitzHugh-Nagumo equations have been used as a caricature of the Hodgkin-Huxley equations of neuron firing and to capture, qualitatively, the general properties of an excitable membrane. In this paper, we utilize a modified version of the FitzHugh-Nagumo equations to model the spatial propagation of neuron firing; we assume that this propagation is (at least, partially) caused by the cross-diffusion connection between the potential and recovery variables. We show that the cross-diffusion version of the model, be- sides giving rise to the typical fast traveling wave solution exhibited in the original "diffusion" FitzHugh-Nagumo equations, additionally gives rise to a slow traveling wave solution. We analyze all possible traveling wave solutions of the model and show that there exists a threshold of the cross-diffusion coefficient (for a given speed of propagation), which bounds the area where "normal" impulse propagation is possible.

Population models with singular equilibrium.

A class of models of biological population and communities with a singular equilibrium at the origin is analyzed; it is shown that these models can possess a dynamical regime of deterministic extinction, which is crucially important from the biological standpoint. This regime corresponds to the presence of a family of homoclinics to the origin, so-called elliptic sector. The complete analysis of possible topological structures in a neighborhood of the origin, as well as asymptotics to orbits tending to this point, is given. An algorithmic approach to analyze system behavior with parameter changes is presented. The developed methods and algorithm are applied to existing mathematical models of biological systems. In particular, we analyze a model of anticancer treatment with oncolytic viruses, a parasite-host interaction model, and a model of Chagas' disease.

Mathematical modeling of tumor therapy with oncolytic viruses: regimes with complete tumor elimination within the framework of deterministic models.

BACKGROUND: Oncolytic viruses that specifically target tumor cells are promising anti-cancer therapeutic agents. The interaction between an oncolytic virus and tumor cells is amenable to mathematical modeling using adaptations of techniques employed previously for modeling other types of virus-cell interaction. RESULTS: A complete parametric analysis of dynamic regimes of a conceptual model of anti-tumor virus therapy is presented. The role and limitations of mass-action kinetics are discussed. A functional response, which is a function of the ratio of uninfected to infected tumor cells, is proposed to describe the spread of the virus infection in the tumor. One of the main mathematical features of ratio-dependent models is that the origin is a complicated equilibrium point whose characteristics determine the main properties of the model. It is shown that, in a certain area of parameter values, the trajectories of the model form a family of homoclinics to the origin (so-called elliptic sector). Biologically, this means that both infected and uninfected tumor cells can be eliminated with time, and complete recovery is possible as a result of the virus therapy within the framework of deterministic models. CONCLUSION: Our model, in contrast to the previously published models of oncolytic virus-tumor interaction, exhibits all possible outcomes of oncolytic virus infection, i.e., no effect on the tumor, stabilization or reduction of the tumor load, and complete elimination of the tumor. The parameter values that result in tumor elimination, which is, obviously, the desired outcome, are compatible with some of the available experimental data. REVIEWERS: This article was reviewed by Mikhail Blagosklonny, David Krakauer, Erik Van Nimwegen, and Ned Wingreen. OPEN PEER REVIEW: Reviewed by Mikhail Blagosklonny, David Krakauer, Erik Van Nimwegen, and Ned Wingreen. For the full reviews, please go to the Reviewers' comments section.

Modeling genome evolution with a diffusion approximation of a birth-and-death process.

MOTIVATION: In our previous studies, we developed discrete-space birth, death and innovation models (BDIMs) of genome evolution. These models explain the origin of the characteristic Pareto distribution of paralogous gene family sizes in genomes, and model parameters that provide for the evolution of these distributions within a realistic time frame have been identified. However, extracting the temporal dynamics of genome evolution from discrete-space BDIM was not technically feasible. We were interested in obtaining dynamic portraits of the genome evolution process by developing a diffusion approximation of BDIM. RESULTS: The diffusion version of BDIM belongs to a class of continuous-state models whose dynamics is described by the Fokker-Plank equation and the stationary solution could be any specified Pareto function. The diffusion models have time-dependent solutions of a special kind, namely, generalized self-similar solutions, which describe the transition from one stationary distribution of the system to another; this provides for the possibility of examining the temporal dynamics of genome evolution. Analysis of the generalized self-similar solutions of the diffusion BDIM reveals a biphasic curve of genome growth in which the initial, relatively short, self-accelerating phase is followed by a prolonged phase of slow deceleration. This evolutionary dynamics was observed both when genome growth started from zero and proceeded via innovation (a potential model of primordial evolution), and when evolution proceeded from one stationary state to another. In biological terms, this regime of evolution can be tentatively interpreted as a punctuated-equilibrium-like phenomenon whereby evolutionary transitions are accompanied by rapid gene amplification and innovation, followed by slow relaxation to a new stationary state.

Gene family evolution: an in-depth theoretical and simulation analysis of non-linear birth-death-innovation models.

BACKGROUND: The size distribution of gene families in a broad range of genomes is well approximated by a generalized Pareto function. Evolution of ensembles of gene families can be described with Birth, Death, and Innovation Models (BDIMs). Analysis of the properties of different versions of BDIMs has the potential of revealing important features of genome evolution. RESULTS: In this work, we extend our previous analysis of stochastic BDIMs. In addition to the previously examined rational BDIMs, we introduce potentially more realistic logistic BDIMs, in which birth/death rates are limited for the largest families, and show that their properties are similar to those of models that include no such limitation. We show that the mean time required for the formation of the largest gene families detected in eukaryotic genomes is limited by the mean number of duplications per gene and does not increase indefinitely with the model degree. Instead, this time reaches a minimum value, which corresponds to a non-linear rational BDIM with the degree of approximately 2.7. Even for this BDIM, the mean time of the largest family formation is orders of magnitude greater than any realistic estimates based on the timescale of life's evolution. We employed the embedding chains technique to estimate the expected number of elementary evolutionary events (gene duplications and deletions) preceding the formation of gene families of the observed size and found that the mean number of events exceeds the family size by orders of magnitude, suggesting a highly dynamic process of genome evolution. The variance of the time required for the formation of the largest families was found to be extremely large, with the coefficient of variation >> 1. This indicates that some gene families might grow much faster than the mean rate such that the minimal time required for family formation is more relevant for a realistic representation of genome evolution than the mean time. We determined this minimal time using Monte Carlo simulations of family growth from an ensemble of simultaneously evolving singletons. In these simulations, the time elapsed before the formation of the largest family was much shorter than the estimated mean time and was compatible with the timescale of evolution of eukaryotes. CONCLUSIONS: The analysis of stochastic BDIMs presented here shows that non-linear versions of such models can well approximate not only the size distribution of gene families but also the dynamics of their formation during genome evolution. The fact that only higher degree BDIMs are compatible with the observed characteristics of genome evolution suggests that the growth of gene families is self-accelerating, which might reflect differential selective pressure acting on different genes.

Birth and death of protein domains: a simple model of evolution explains power law behavior.

BACKGROUND: Power distributions appear in numerous biological, physical and other contexts, which appear to be fundamentally different. In biology, power laws have been claimed to describe the distributions of the connections of enzymes and metabolites in metabolic networks, the number of interactions partners of a given protein, the number of members in paralogous families, and other quantities. In network analysis, power laws imply evolution of the network with preferential attachment, i.e. a greater likelihood of nodes being added to pre-existing hubs. Exploration of different types of evolutionary models in an attempt to determine which of them lead to power law distributions has the potential of revealing non-trivial aspects of genome evolution. RESULTS: A simple model of evolution of the domain composition of proteomes was developed, with the following elementary processes: i) domain birth (duplication with divergence), ii) death (inactivation and/or deletion), and iii) innovation (emergence from non-coding or non-globular sequences or acquisition via horizontal gene transfer). This formalism can be described as a birth, death and innovation model (BDIM). The formulas for equilibrium frequencies of domain families of different size and the total number of families at equilibrium are derived for a general BDIM. All asymptotics of equilibrium frequencies of domain families possible for the given type of models are found and their appearance depending on model parameters is investigated. It is proved that the power law asymptotics appears if, and only if, the model is balanced, i.e. domain duplication and deletion rates are asymptotically equal up to the second order. It is further proved that any power asymptotic with the degree not equal to -1 can appear only if the hypothesis of independence of the duplication/deletion rates on the size of a domain family is rejected. Specific cases of BDIMs, namely simple, linear, polynomial and rational models, are considered in details and the distributions of the equilibrium frequencies of domain families of different size are determined for each case. We apply the BDIM formalism to the analysis of the domain family size distributions in prokaryotic and eukaryotic proteomes and show an excellent fit between these empirical data and a particular form of the model, the second-order balanced linear BDIM. Calculation of the parameters of these models suggests surprisingly high innovation rates, comparable to the total domain birth (duplication) and elimination rates, particularly for prokaryotic genomes. CONCLUSIONS: We show that a straightforward model of genome evolution, which does not explicitly include selection, is sufficient to explain the observed distributions of domain family sizes, in which power laws appear as asymptotic. However, for the model to be compatible with the data, there has to be a precise balance between domain birth, death and innovation rates, and this is likely to be maintained by selection. The developed approach is oriented at a mathematical description of evolution of domain composition of proteomes, but a simple reformulation could be applied to models of other evolving networks with preferential attachment.