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.

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.

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.

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.

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.

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.

Simple stochastic birth and death models of genome evolution: was there enough time for us to evolve?

MOTIVATION: The distributions of many genome-associated quantities, including the membership of paralogous gene families can be approximated with power laws. We are interested in developing mathematical models of genome evolution that adequately account for the shape of these distributions and describe the evolutionary dynamics of their formation. RESULTS: We show that simple stochastic models of genome evolution lead to power-law asymptotics of protein domain family size distribution. These models, called Birth, Death and Innovation Models (BDIM), represent a special class of balanced birth-and-death processes, in which domain duplication and deletion rates are asymptotically equal up to the second order. The simplest, linear BDIM shows an excellent fit to the observed distributions of domain family size in diverse prokaryotic and eukaryotic genomes. However, the stochastic version of the linear BDIM explored here predicts that the actual size of large paralogous families is reached on an unrealistically long timescale. We show that introduction of non-linearity, which might be interpreted as interaction of a particular order between individual family members, allows the model to achieve genome evolution rates that are much better compatible with the current estimates of the rates of individual duplication/loss events.

The structure of the protein universe and genome evolution.

Despite the practically unlimited number of possible protein sequences, the number of basic shapes in which proteins fold seems not only to be finite, but also to be relatively small, with probably no more than 10,000 folds in existence. Moreover, the distribution of proteins among these folds is highly non-homogeneous -- some folds and superfamilies are extremely abundant, but most are rare. Protein folds and families encoded in diverse genomes show similar size distributions with notable mathematical properties, which also extend to the number of connections between domains in multidomain proteins. All these distributions follow asymptotic power laws, such as have been identified in a wide variety of biological and physical systems, and which are typically associated with scale-free networks. These findings suggest that genome evolution is driven by extremely general mechanisms based on the preferential attachment principle.

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.