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.

How trait distributions evolve in populations with parametric heterogeneity.

We consider the problem of determining the time evolution of a trait distribution in a mathematical model of non-uniform populations with parametric heterogeneity. This means that we consider only heterogeneous populations in which heterogeneity is described by an individual specific parameter that differs in general from individual to individual, but does not change with time for the whole lifespan of this individual. Such a restriction allows obtaining a number of simple and yet important analytical results. In particular we show that initial assumptions on time-dependent behavior of various characteristics, such as the mean, variance, or coefficient of variation, restrict severely possible choices for the exact form of the trait distribution. This fact must be taken into account for both model formulation and, especially, for testing theoretical models against available real world data. We illustrate our findings by in-depth analysis of the variance evolution and specific examples from population ecology and mathematical epidemiology. We also reanalyze a well known mathematical model for gypsy moth population and show that the knowledge of how trait distributions evolve allows producing oscillatory behaviors for highly heterogeneous populations.

Natural Selection Between Two Games with Applications to Game Theoretical Models of Cancer.

Evolutionary game theory has been used extensively to study single games as applied to cancer, including in the context of metabolism, development of resistance, and even games between tumor and treatment. However, the situation when several games are being played against each other at the same time has not yet been investigated. Here, we describe a mathematical framework for analyzing natural selection not just between strategies, but between games. We provide theoretical analysis of situations of natural selection between the games of Prisoner's dilemma and Hawk-Dove, and demonstrate that while the dynamics of cooperators and defectors within their respective games is as expected, the distribution of games changes over time due to natural selection. We also investigate the question of mutual invasibility of games with respect to different strategies and different initial population composition. We conclude with a discussion of how the proposed approach can be applied to other games in cancer, such as motility versus stability strategies that underlie the process of metastatic invasion.

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.

From Experiment to Theory: What Can We Learn from Growth Curves?

Finding an appropriate functional form to describe population growth based on key properties of a described system allows making justified predictions about future population development. This information can be of vital importance in all areas of research, ranging from cell growth to global demography. Here, we use this connection between theory and observation to pose the following question: what can we infer about intrinsic properties of a population (i.e., degree of heterogeneity, or dependence on external resources) based on which growth function best fits its growth dynamics? We investigate several nonstandard classes of multi-phase growth curves that capture different stages of population growth; these models include hyperbolic-exponential, exponential-linear, exponential-linear-saturation growth patterns. The constructed models account explicitly for the process of natural selection within inhomogeneous populations. Based on the underlying hypotheses for each of the models, we identify whether the population that it best fits by a particular curve is more likely to be homogeneous or heterogeneous, grow in a density-dependent or frequency-dependent manner, and whether it depends on external resources during any or all stages of its development. We apply these predictions to cancer cell growth and demographic data obtained from the literature. Our theory, if confirmed, can provide an additional biomarker and a predictive tool to complement experimental research.

Non-linearity and heterogeneity in modeling of population dynamics.

The study of population growth reveals that the behaviors that follow the power law appear in numerous biological, demographical, ecological, physical and other contexts. Parabolic models appear to be realistic approximations of real-life replicator systems, while hyperbolic models were successfully applied to problems of global demography and appear relevant in quasispecies and hypercycle modeling. Nevertheless, it is not always clear why non-exponential growth is observed empirically and what possible origins of the non-exponential models are. In this paper the power equation is considered within the frameworks of inhomogeneous population models; it is proven that any power equation describes the total population size of a frequency-dependent model with Gamma-distributed Malthusian parameter. Additionally, any super-exponential equation describes the dynamics of inhomogeneous Malthusian density-dependent population model. All statistical characteristics of the underlying inhomogeneous models are computed explicitly. The results of this analysis show that population heterogeneity can be a reasonable explanation for power law accurately describing total population growth.

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.

Mixed strategies and natural selection in resource allocation.

An appropriate choice of strategy for resource allocation may frequently determine whether a population will be able to survive under the conditions of severe resource limitations. Here we focus on two classes of strategies allocation of resources towards rapid proliferation, or towards slower proliferation but increased physiological and environmental maintenance. We propose a generalized framework, where individuals within a population can use either strategy in different proportion for utilization of a common dynamical resource in order to maximize their fitness. We use the model to address two major questions, namely, whether either strategy is more likely to be selected for as a result of natural selection, and, if one allows for the possibility of resource over-consumption, whether either strategy is preferable for avoiding population collapse due to resource exhaustion. Analytical and numerical results suggest that the ultimate choice of strategy is determined primarily by the initial distribution of individuals in the population, and that while investment in physiological and environmental maintenance is a preferable strategy in a homogeneous population, no generalized prediction can be made about heterogeneous populations.

Parabolic replicator dynamics and the principle of minimum Tsallis information gain.

BACKGROUND: Non-linear, parabolic (sub-exponential) and hyperbolic (super-exponential) models of prebiological evolution of molecular replicators have been proposed and extensively studied. The parabolic models appear to be the most realistic approximations of real-life replicator systems due primarily to product inhibition. Unlike the more traditional exponential models, the distribution of individual frequencies in an evolving parabolic population is not described by the Maximum Entropy (MaxEnt) Principle in its traditional form, whereby the distribution with the maximum Shannon entropy is chosen among all the distributions that are possible under the given constraints. We sought to identify a more general form of the MaxEnt principle that would be applicable to parabolic growth. RESULTS: We consider a model of a population that reproduces according to the parabolic growth law and show that the frequencies of individuals in the population minimize the Tsallis relative entropy (non-additive information gain) at each time moment. Next, we consider a model of a parabolically growing population that maintains a constant total size and provide an "implicit" solution for this system. We show that in this case, the frequencies of the individuals in the population also minimize the Tsallis information gain at each moment of the 'internal time" of the population. CONCLUSIONS: The results of this analysis show that the general MaxEnt principle is the underlying law for the evolution of a broad class of replicator systems including not only exponential but also parabolic and hyperbolic systems. The choice of the appropriate entropy (information) function depends on the growth dynamics of a particular class of systems. The Tsallis entropy is non-additive for independent subsystems, i.e. the information on the subsystems is insufficient to describe the system as a whole. In the context of prebiotic evolution, this "non-reductionist" nature of parabolic replicator systems might reflect the importance of group selection and competition between ensembles of cooperating replicators.

Preventing the tragedy of the commons through punishment of over-consumers and encouragement of under-consumers.

The conditions that can lead to the exploitative depletion of a shared resource, i.e., the tragedy of the commons, can be reformulated as a game of prisoner's dilemma: while preserving the common resource is in the best interest of the group, over-consumption is in the interest of each particular individual at any given point in time. One way to try and prevent the tragedy of the commons is through infliction of punishment for over-consumption and/or encouraging under-consumption, thus selecting against over-consumers. Here, the effectiveness of various punishment functions in an evolving consumer-resource system is evaluated within a framework of a parametrically heterogeneous system of ordinary differential equations (ODEs). Conditions leading to the possibility of sustainable coexistence with the common resource for a subset of cases are identified analytically using adaptive dynamics; the effects of punishment on heterogeneous populations with different initial composition are evaluated using the reduction theorem for replicator equations. Obtained results suggest that one cannot prevent the tragedy of the commons through rewarding of under-consumers alone--there must also be an implementation of some degree of punishment that increases in a nonlinear fashion with respect to over-consumption and which may vary depending on the initial distribution of clones in the population.

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.

The universal distribution of evolutionary rates of genes and distinct characteristics of eukaryotic genes of different apparent ages.

The evolutionary rates of protein-coding genes in an organism span, approximately, 3 orders of magnitude and show a universal, approximately log-normal distribution in a broad variety of species from prokaryotes to mammals. This universal distribution implies a steady-state process, with identical distributions of evolutionary rates among genes that are gained and genes that are lost. A mathematical model of such process is developed under the single assumption of the constancy of the distributions of the propensities for gene loss (PGL). This model predicts that genes of different ages, that is, genes with homologs detectable at different phylogenetic depths, substantially differ in those variables that correlate with PGL. We computationally partition protein-coding genes from humans, flies, and Aspergillus fungus into age classes, and show that genes of different ages retain the universal log-normal distribution of evolutionary rates, with a shift toward higher rates in "younger" classes but also with a substantial overlap. The only exception involves human primate-specific genes that show a heavy tail of rapidly evolving genes, probably owing to gene annotation artifacts. As predicted, the gene age classes differ in characteristics correlated with PGL. Compared with "young" genes (e.g., mammal-specific human ones), "old" genes (e.g., eukaryote-specific), on average, are longer, are expressed at a higher level, possess a higher intron density, evolve slower on the short time scale, and are subject to stronger purifying selection. Thus, genome evolution fits a simple model with approximately uniform rates of gene gain and loss, without major bursts of genomic innovation.

"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: effects of parametric heterogeneity on cell dynamics.

BACKGROUND: One of the mechanisms that ensure cancer robustness is tumor heterogeneity, and its effects on tumor cells dynamics have to be taken into account when studying cancer progression. There is no unifying theoretical framework in mathematical modeling of carcinogenesis that would account for parametric heterogeneity. RESULTS: Here we formulate a modeling approach that naturally takes stock of inherent cancer cell heterogeneity and illustrate it with a model of interaction between a tumor and an oncolytic virus. We show that several phenomena that are absent in homogeneous models, such as cancer recurrence, tumor dormancy, and others, appear in heterogeneous setting. We also demonstrate that, within the applied modeling framework, to overcome the adverse effect of tumor cell heterogeneity on the outcome of cancer treatment, a heterogeneous population of an oncolytic virus must be used. Heterogeneity in parameters of the model, such as tumor cell susceptibility to virus infection and the ability of an oncolytic virus to infect tumor cells, can lead to complex, irregular evolution of the tumor. Thus, quasi-chaotic behavior of the tumor-virus system can be caused not only by random perturbations but also by the heterogeneity of the tumor and the virus. CONCLUSION: The modeling approach described here reveals the importance of tumor cell and virus heterogeneity for the outcome of cancer therapy. It should be straightforward to apply these techniques to mathematical modeling of other types of anticancer therapy. REVIEWERS: Leonid Hanin (nominated by Arcady Mushegian), Natalia Komarova (nominated by Orly Alter), and David Krakauer.

Biological applications of the theory of birth-and-death processes.

In this review, we discuss applications of the theory of birth-and-death processes to problems in biology, primarily, those of evolutionary genomics. The mathematical principles of the theory of these processes are briefly described. Birth-and-death processes, with some straightforward additions such as innovation, are a simple, natural and formal framework for modeling a vast variety of biological processes such as population dynamics, speciation, genome evolution, including growth of paralogous gene families and horizontal gene transfer and somatic evolution of cancers. We further describe how empirical data, e.g. distributions of paralogous gene family size, can be used to choose the model that best reflects the actual course of evolution among different versions of birth-death-and-innovation models. We conclude that birth-and-death processes, thanks to their mathematical transparency, flexibility and relevance to fundamental biological processes, are going to be an indispensable mathematical tool for the burgeoning field of systems biology.

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.

Mathematical modeling of evolution of horizontally transferred genes.

We describe a stochastic birth-and-death model of evolution of horizontally transferred genes in microbial populations. The model is a generalization of the stochastic model described by Berg and Kurland and includes five parameters: the rate of mutational inactivation, selection coefficient, invasion rate (i.e., rate of arrival of a novel sequence from outside of the recipient population), within-population horizontal transmission ("infection") rate, and population size. The model of Berg and Kurland included four parameters, namely, mutational inactivation, selection coefficient, population size, and "infection." However, the effect of "infection" was disregarded in the interpretation of the results, and the overall conclusion was that horizontally acquired sequences can be fixed in a population only when they confer a substantial selective advantage onto the recipient and therefore are subject to strong positive selection. Analysis of the present model in different domains of parameter values shows that, as long as the rate of within-population horizontal transmission is comparable to the mutational inactivation rate and there is even a low rate of invasion, horizontally acquired sequences can be fixed in the population or at least persist for a long time in a substantial fraction of individuals in the population even when they are neutral or slightly deleterious. The available biological data strongly suggest that intense within-population and even between-populations gene flows are realistic for at least some prokaryotic species and environments. Therefore, our modeling results are compatible with the notion of a pivotal role of horizontal gene transfer in the evolution of prokaryotes.

A simple epidemic model with surprising dynamics.

A simple model incorporating demographic and epidemiological processes is explored. Four re-parameterized quantities the basic demographic reproductive number (R(d)), the basic epidemiological reproductive number (R(0)), the ratio (v) between the average life spans of susceptible and infective class, and the relative fecundity of infectives (theta), are utilized in qualitative analysis. Mathematically, non-analytic vector fields are handled by blow-up transformations to carry out a complete and global dynamical analysis. A family of homoclinics is found, suggesting that a disease outbreak would be ignited by a tiny number of infectious individuals.