 Research
 Open Access
 Published:
Pseudochaotic oscillations in CRISPRvirus coevolution predicted by bifurcation analysis
Biology Direct volume 9, Article number: 13 (2014)
Abstract
Background
The CRISPRCas 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 CRISPRCas has sharply increased the interest in the diversity and evolution of these systems. Comparative genomic data indicate that during evolution of prokaryotes CRISPRCas loci are lost and acquired via horizontal gene transfer at high rates. Mathematical modeling and initial experimental studies of CRISPRcarrying microbes and viruses reveal complex coevolutionary dynamics.
Results
We performed a bifurcation analysis of models of coevolution of viruses and microbial host that possess CRISPRCas hereditary adaptive immunity systems. The analyzed Malthusian and logistic models display complex, and in particular, quasichaotic oscillation regimes that have not been previously observed experimentally or in agentbased models of the CRISPRmediated immunity. The key factors for the appearance of the quasichaotic oscillations are the nonlinear 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 CRISPRhost coevolution model predicts complex regimes including quasichaotic oscillations. The quasichaotic regimes of virushost coevolution are likely to be biologically relevant given the evolutionary instability of the CRISPRCas loci revealed by comparative genomics. The results of this analysis might have implications beyond the CRISPRCas 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.
Background
The arms races between microbes and viruses preying on them often display rich, complex population dynamics [1]. In principle, the dynamics of virusmicrobe interactions is analogous to the classical predator–prey models [2–5] but both microbes and viruses evolve much faster than animals such that virushost interactions change on a scale that may be amenable to direct laboratory study. One of the adaptation mechanisms employed by hosts to curb viruses is the CRISPRCas (Clustered Regularly Interspaced Short Palindromic RepeatsCRISPR associated proteins), a recently discovered adaptive immunity system that is present in the great majority of Archaea and many bacteria [6–12]. Microbes create heritable memory of viruses that attack them by inserting virusderived spacers into CRISPR repeat cassettes, thus following the Lamarckian modality of evolution that dramatically accelerates adaptation [13]. The rapid adaptation through the activity of CRISPRCas is possible because this system engenders heritable genetic changes that are directly beneficial for the archaeon or bacterium in the face of a specific environmental challenge (a virus), in contrast to the random, undirected mutations in the Darwinian evolutionary framework [14]. The CRISPRCas systems are increasingly used as powerful, versatile tools for genomic engineering tools which sharply increases the interest in their diversity and evolution [15–20].
General considerations suggest that the population dynamics of virushost coevolution should be dominated by periodic selective sweeps alternating between the host (when it "discovers" a resistance mutation or acquires immunity against the dominant virus lineage) and the virus (when an immunity escape mutation occurs), similar to the case of rapidly evolving human viruses [21]. Indeed, such behavior has been observed in simple LotkaVolterra type models of phagebacteria coevolution [22, 23] as well as in direct evolutionary experiments [23]. However, direct and indirect population studies reveal much more complex behaviors of the actual populations that usually does not involve strain dominance and instead displays longterm persistence of multiple lineages of both the microbial hosts and the viruses [24–27].
Several detailed agentbased models of coevolution between viruses and CRISPRCascarrying hosts have been developed and analyzed [26, 28–35]. The agentbased models incorporate the salient features of the CRISPRCas system such as the existence of the CRISPR cassette with virusderived spacers, immunity conferred to a host by spacers that match the attacking virus and acquisition of new spacers as a result of failed virus infections. These models allow one to reproduce many aspects of the observed behavior of coevolving virushost systems and predict conditions required for the evolutionary maintenance of the CRISPRCas immunity, such as a threshold of viral diversity [28].
Agentbased models provide for the exploration of interactions of arbitrary complexity and naturally incorporate the desired level of granularity (e.g. individualbased or lineagebased models) and the stochasticity of the processes involved. However, such models typically possess a highdimensional parameter space that cannot be explored in full, so that not all potential regimes, some of which could be biologically relevant, are captured. In contrast, mathematical models based on systems of differential equations are limited in complexity and are inherently less realistic (at the very least because they approximate reality with infinitely small deterministic changes) but when analytically tractable, permit a full and rigorous analysis of all possible behaviors.
Here we describe LotkaVolterra type models of interaction between a host with a heritable adaptive immunity system, such as CRISPRCas, and a virus that escapes the immunity via implicit accumulation of mutations which is implemented as gradual immunity decay. We construct "minimal" analytical models which capture qualitatively the basic regimes of the CRISPRCas system behaviors previously found experimentally and through the agentbased modeling. We explore the full spectrum of possible behaviors of this virushost system, compare the results with those of a more detailed agentbased model [32], and describe a previously unnoticed regime of quasichaotic oscillations.
Results
Threecomponent CRISPR population dynamics: Malthusian and logistic versions
In order to construct "minimal" analytical models that would qualitatively capture the basic regimes of the CRISPRCas system behaviors we analyzed, as a preliminary step, a twocomponent Volterratype model which describes the dynamics of virushost system and consider immunity static within the timescale of the model. Our analysis shows that twocomponent models display simple dynamical regimes (see Methods for details). However, the CRISPRCas system of adaptive immunity, which this work seeks to model, cannot be expected to follow these simple regimes given its dynamic evolution that involves rapid acquisition of immunity to a particular virus or plasmid along with loss and gain of entire CRISPRCas loci [36]. Thus, more realistic models should take into account that adaptive immunity systems are characterized by the existence of immune memory that enhances the response in individuals encountering a familiar challenge [37, 38]. Originally nonimmune ("native") individuals can acquire the adaptive response capacity that persists at timescales relevant for our model. Thus, we introduce two categories of hosts, immune and nonimmune, that interact differently with the virus and exchange individuals via immunity acquisition and decay.
Let x(t) be the density of immune hosts with the immunity p with respect to viruses, 0 ≤ p ≤ 1; y(t) be the density of sensitive hosts with immunity s, 0 ≤ s ≤ p ≤ 1; z(t) be the density of viruses. We consider the following 3component model:
Here l is the immunity decay rate (x → y flow), e is the immunity acquisition rate; d is the death rate of viruses, M is the virus reproduction rate; b is the encounter rate coefficient; the growth rates of immune and sensitive hosts are equal to 1.
Model (1) for a = 0 describes a situation when immune and sensitive hosts in the absence of viruses grow according to the Malthusian model. If a > 0, the model (1) describes a more realistic situation when both classes of hosts grow according to the logistics model.
Coordinates of equilibria of (1) solve the system:
Model (1) has equilibrium O(0, 0, 0) in all cases. The logistic model with a > 0 has one more equilibrium A(0,1/a,0), which corresponds to existence of only nonimmune hosts.
In Methods the following assertions are proven.
Statement 1.

(1)
Equilibrium O(x = y = z = 0) is unstable for all parameter values;

(2)
If a > 0 then equilibrium A(0,1/a,0) is stable for $\mathit{M}<\frac{\mathit{\text{da}}+\mathit{\text{bs}}}{\mathit{b}\left(1\mathit{s}\right)}$ and unstable for $\mathit{M}>\frac{\mathit{\text{da}}+\mathit{\text{bs}}}{\mathit{b}\left(1\mathit{s}\right)}.$
If the immunity p in model (1) is a constant, then, additionally, the model may have either nontrivial stable equilibrium, or may demonstrate periodic oscillations. Detail description of possible behaviors of the model with constant p is given in Methods.
More realistically, the immunity p is not a constant but depends on the density of the virus, p = p(z) (0 ≤ p ≤ 1). Below we consider the latter case, in agreement with empirical observations and computer simulations. Specifically, it has been shown that CRISPRCas systems are (nearly) ubiquitous in archaeal and bacterial hyperthermophiles are present in less than half of the available mesophile genomes [28, 32, 39, 40]. Analysis of agentbased models of virushost coevolution suggest that this distinction stems from the fact that hyperthermophiles face lower virus loads and diversity than mesophiles providing for higher efficacy of CRISPRCas [28].
Our aim is to find all stable modes of the model at different values of the model parameters and to describe the transitions from one mode to another when parameters vary; by other words, we want to construct the bifurcation diagram of the model. It is natural to suppose that p = p(z) monotonically decreases and tends to the immunity s of sensitive hosts at large z. From now on we consider:
where k, s are constants, 0 < s < 1, k > 0. Under equation (3), immunity is a monotonically declining function of the virus amount that tends to a constant, maximum p (maximally efficient adaptive immunity) when z tends to zero, and tends to s (no adaptive immunity, innate immunity only) when z tends to infinity.
Let us consider model (1) with the immunity p = p(z) defined by (3). We do not attempt a complete analysis of this model but rather seek to identify stable modes and most interesting dynamical behaviors.
We start with the Malthusian version when a = 0. It has nontrivial equilibrium B_{ e }(x_{ e }, y_{ e }, z_{ e }) such that x = x_{ e }, y = y_{ e } are expressed via z = z_{ e }:
where z = z_{ e } solves the equation
In particular, we show that for a wide domain of the parameter values, the model demonstrates nonperiodic oscillation of all three variables. The following assertions are valid (see Methods).
Statement 2. For a wide range of (fixed) parameters l, e, s, b, 0 < k ≤ 1, system (1), (3) with a = 0 has only one positive and unstable equilibrium B_{ e }(x_{ e }, y_{ e }, z_{ e }) under condition $\mathit{p}({\mathit{z}}_{\mathit{e}})<\frac{\mathit{M}}{\mathit{M}+1};$the coordinates (x_{ e }, y_{ e }, z_{ e }) of this equilibrium satisfy (4), (5).
Trajectories of the model show quasichaotic behavior for a broad range of the model parameters. Consider in detail the behavior of the model solutions depending on the values of the parameters l, e, and M. For l > e, typical trajectories starting close to the equilibrium point are shown at Figure 1. Initially, all variables show almost periodic oscillations with increasing amplitudes. Then, as the amplitudes become large enough, the behavior of the trajectories changes sharply and the oscillations become (quasi)chaotic; if the initial point is far from the equilibrium, then the (quasi)chaotic oscillations are observed from the very beginning. When l < e the behavior of the trajectories is similar. The difference is that the fraction of immune hosts, x, in case l < e is greater than it is in case l > e. Again, if the initial point is taken far from the equilibrium, then the (quasi)chaotic oscillations are observed from the very beginning, similar to Figure 1. These types of behavior are observed in a wide area of values of the parameter M, 1 < M < 1000. Notice, that when M increases, the maximum values of x,y decrease whereas z does not depend on M. This effect does not seem to have a plausible biological interpretation (viruses cannot exist if the hosts go extinct), indicative of apparent limitations of the model.
Let us consider a more realistic version of 3component system (1) with logistic growth of hosts and the immunity p(z) given by (3). Again, we do not attempt a complete analysis of this model but rather seek to identify stable modes and most interesting dynamical behaviors and to compare them to those of the Malthusian version of the model. In particular, we show that the model displays nontrivial stable equilibria, stable oscillations, and quasichaotic oscillations of all three variables. More precisely, there exists a critical value of the virus birth rate M^{cr} at fixed values of all other parameters such that the system tends to a stable equilibrium when M < M^{cr}, but if M > M^{cr}, then small periodic oscillations appear due to the Hopf bifurcation and then, if M > > M^{cr}, the system transits to a regime of (quasi)chaotic oscillations.
Coordinates of any nontrivial equilibrium B(x_{ e }, y_{ e }, z_{ e }) should satisfy system (2); taking a = 1 for simplicity, we present system (2) in the form:
It follows from the first two equations of (6) that
The solution to the last equation is $\phantom{\rule{0.25em}{0ex}}\mathit{x}+\mathit{y}=\frac{1+\mathit{M}\mathit{\text{bz}}+\sqrt{{\left(1+\mathit{\text{bz}}\mathit{M}\right)}^{2}4\mathit{d}\left(1+\mathit{M}\right)\mathit{z}}}{2\left(1+\mathit{M}\right)}\le \frac{1+\mathit{M}\mathit{\text{bz}}}{1+\mathit{M}}.$ So,
It means that all possible nontrivial equilibria of the model are placed in a bounded area of the variable values; the amount of hosts is comparatively small (<1) and the amount of viruses is restricted by the virus reproduction rate M.
Solving system (6) we find coordinates of nontrivial equilibrium B(x_{ e }, y_{ e }, z_{ e }) such that x = x_{ e }, y = y_{ e } are expressed via z = z_{ e }
and z coordinate solves the equation:
Analysis of the system (1) with a > 0, (3) showed that there exists such threshold M^{cr} (depending on the model parameters) that the equilibrium B is stable if M < M^{cr}; when M increases and intersects the threshold M = M^{cr}, the equilibrium B loses stability and a stable limit cycle appears in the system. We summarize the results of this analysis in
Statement 3. For a wide range of (fixed) parameters l, e, s, b, 0 < k ≤ 1, model (1), (3) with a = 1 has a stable equilibrium for 0 < M < M^{cr}and stable oscillations for M > M^{cr}, which appear due to supercritical Hopf bifurcation at M = M^{cr}(l, e, s, b).
See Methods for sketch of the proof.
Examples of the evolution of trajectories and phase portraits as the value of M increases are given in Figures 2, 3, 4, 5. Figure 2 shows trajectories of the system that tend to a stable equilibrium when the parameter M is below the bifurcation threshold. Figure 3 demonstrates that the system arrives at a stable limit cycle when M intersects the bifurcation threshold, and Figure 4 shows that the established regime at the steady state is very close to periodic oscillations.
Oscillations that are observed in Figure 4, with M above the bifurcation threshold M^{cr}, are very close to but not perfectly periodic. This peculiarity of the trajectories can be explained in the following way. The Hopf theorem for a 3dimension system states (see, e.g., [41], ch.5) that there exists such onetoone transformation of the initial variables x → x*, y → y*, z → z* that two of new variables (say, x* and y*) show stable periodic oscillations but the third one does not and is governed by a separate equation. For this reason, the trajectories of initial variables, which are functions of x*, y*, z*, may be nonperiodic and may show not exactly periodic and even quasi chaotic oscillations.
It should be emphasized that, when the value of parameter M crosses the bifurcation boundary, the qualitative behavior of the system changes sharply: as M increases, the behavior of the system becomes more and more "chaotic", with increasing nonregularity of the shapes of the trajectories (Figure 5). At large M, phase curves fill a surface in the (x, y, z) space (see Figure 6, left panel, in contrast to Figure 4), and the phase portrait of the system reveals sharp, quasichaotic oscillations (see Figure 6, right panel).
Discussion
The bifurcation analysis of the virushost system described here reveals complex, and in particular quasichaotic, oscillation regimes that so far have not been previously observed experimentally or in agentbased models of the CRISPRmediated adaptive immunity [26, 28–35]. The patchy distribution of the CRISPRCas systems across the microbial diversity, and even within relatively narrow groups of bacteria [10, 42, 43], implies complex dynamics of virushost coevolution. Here we show that, in order to detect such nontrivial coevolutionary regimes, the model has to be sufficiently complex, or put another way, less unrealistic than toy models that deal with a single host population and a simple dependence (or independence) of immunity on the amount of the virus, such as the twocomponent models examined here. The key factors for the appearance of quasichaotic oscillations are the nonlinear dependence of the immunity on the amount of viruses and the threedimensional phase space of the model which divides the hosts into the immune and susceptible populations, so that the system consists of three components (Table 1). In a conceptually similar manner, chaotic behavior has been observed in a model of bacteriophagehost interaction where the complexity of the system is introduced through inclusion of the timedelayed formation of consumable bacterial debris as a result of cell lysis [44, 45].
It should be emphasized that the threecomponent models require the immune host to persist in the population, hence the results obtained here are only applicable to adaptive immunity that involves stable inheritance, i.e. the Lamarckian mode of evolution [13]. The CRISPRCas systems that function by introducing unique, directed modifications into the host genomes present the most straightforward case of such Lamarckian immunity [13, 34]. Nevertheless, other adaptive immunity systems, in particular the piRNA mechanism of transposon restriction in the animal germline, function on similar principles [46]. Moreover, the siRNA branch of eukaryotic RNA interference, which is nearly universal among eukaryotes, also encompasses a substantial heritable component, albeit in this case via the epigenetic route [37, 47–49]. The quasichaotic regime of virushost coevolution described here might be relevant for these immunity systems as well.
Both the Malthusian and the logistic versions of the threecomponent model demonstrate quasichaotic oscillations but they appear via distinct mechanisms (Table 1). The Malthusian version has no stable modes, with unstable equilibria. For a wide area of the parameter values, the trajectories lie in a bounded domain, and for this reason, these trajectories demonstrate (quasi)chaotic behavior for all reasonable values of the virus birth rate M. In contrast, the logistic model has stable modes. There exists a critical value, namely a threshold of the virus birth rate M^{cr} at fixed values of all other parameters, such that the system tends to a stable nontrivial equilibrium at M < M^{cr}. When M > M^{cr}, the qualitative behavior of the system changes sharply. A limit cycle and (almost) periodic oscillations appear. With the further increase of M, the limit cycle turns into a "surface", and the behavior of the system becomes more and more "chaotic", with increasing nonregularity of the shapes of the trajectories. When M > > M^{cr}, the system transits to a regime of quasichaotic oscillations.
Thus, the models described here make concrete, quantitative predictions that can be directly tested using an experimental setup for virushost coevolution [27, 50].
Conclusions
Comparative genomics as well as direct experiments reveal notable evolutionary instability of the CRISPRCas systems [36, 50–52]. It is common for closely related isolates of bacteria or archaea to differ with respect to the presence or absence of CRISPRCas, indicating that this immune system is easily lost and gained via HGT. Rearrangements of the CRISPRCas loci are also extremely widespread among microbes [43, 53]. Given this evolutionary plasticity, complex and in particular quasichaotic regimes of virushost coevolution revealed here appear to be plausible and potentially important for the evolutionary outcomes. In light of the current, rapidly increasing interest in CRISPR research, experimental validation of such regimes could be a realistic prospect.
Methods
Twocomponent model with p = const
The hostvirus dynamics may be considered within the framework of the Volterratype models
If the parameter p is constant and a = 0 then model (M1) is Hamiltonian (conservative) with the Hamiltonian G(x, z) = ln z + d ln x  b(M(1  p)  p)x  b(1  p)z. If $0<\mathit{p}<\frac{\mathit{M}}{\mathit{M}+1}$ then the system has a saddle in the origin O and a center in the equilibrium
If a > 0 (logistic case), then for constant $0<\mathit{p}<\frac{\mathit{M}}{\mathit{M}+1}$ system (M1) has a saddle in the origin, equilibria A(1/a, 0) and $\mathit{B}\left(\mathit{x}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\left(1\mathit{p}\right)\mathit{p}\right)},\mathit{z}=\frac{\mathit{b}\left(\mathit{M}\left(1\mathit{p}\right)\mathit{p}\right)\mathit{\text{ad}}}{{\mathit{b}}^{2}\left(1\mathit{p}\right)\left(\mathit{M}\left(1\mathit{p}\right)\mathit{p}\right)}\right);$ equilibrium B belongs to the positive quadrant (x, z) if b(M(1  p)  p)  ad > 0. In this case, B is a stable node/spiral and A is a saddle , A is a stable node if B is not positive (see, e.g., [54]).
Twocomponent model with p = p(z)
Malthusian version of the model (M1) (a = 0).
This system has equilibrium in the origin (0,0), which is a saddle; z coordinate of any other equilibrium $\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)$ has to be a root of the equation 1  bz(1  p(z)) = 0 and $\overline{\mathit{x}}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\mathit{p}\left(\overline{\mathit{z}}\right)\left(1+\mathit{M}\right)\right)}.$ The point $\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)$ is positive only if $\mathit{p}\left(\overline{\mathit{z}}\right)<\mathit{M}/\left(1+\mathit{M}\right).$
Proposition 1. If p_{ z }(z) < 0, then the Malthusian model has only one nontrivial equilibrium $\mathit{B}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right),$which is an unstable node/spiral for every parameter values.
Proof. The Jacobian of the system in the equilibrium $\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)$ is
and its determinant and trace are:
If the function p(z) decreases monotonically, i.e. p_{ z }(z) < 0, then p(z) and monotonically increasing function $\mathit{h}\left(\mathit{z}\right)=1\frac{1}{\mathit{\text{bz}}}$ can intersect only once. Thus, equation 1  bz(1  p(z)) = 0 has only one root $\overline{\mathit{z}}$, and the Malthusian model has only one nontrivial equilibrium $\mathit{B}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right).$ Next, $\mathit{\text{Det}}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)>0,\mathit{\text{Tr}}\phantom{\rule{0.24em}{0ex}}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)>0$ for positive $\overline{\mathit{x}},\overline{\mathit{z}}$ if p_{ z }(z) < 0. Thus $\mathit{B}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)$ is unstable node or unstable spiral.
Logistic version of model (M1) (a = 1).
Equilibria of system (M1) with a = 1 and p = p(z) defined by (3) are the points (0,0), A(1,0), and $\mathit{B}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)$ where coordinates $\overline{\mathit{x}},\overline{\mathit{z}}$ solve the system
Denote $\mathit{h}\left(\mathit{z}\right)\equiv \frac{\mathit{d}\mathit{\text{bM}}+{\mathit{b}}^{2}\mathit{\text{Mz}}}{\mathit{b}\left(1\mathit{M}+\mathit{\text{bMz}}\right)},$ then $\overline{\mathit{z}}$ is a root of the equation p(z) = h(z). Solutions of this equation are the points of intersection of the curves p(z) and h(z); up to two equilibrium points ${\mathit{B}}_{1}\left({\overline{\mathit{x}}}_{1},{\overline{\mathit{z}}}_{1}\right),{\mathit{B}}_{2}\left({\overline{\mathit{x}}}_{2},{\overline{\mathit{z}}}_{2}\right)$ can appear in the model with parameter variations. Denote J(x, z) the Jacobian of system (M1), (3) with a = 1:
Thus O(0,0) is a saddle and A(1,0) is a stable node for all parameter values.
Now let us consider the system consisting of equation (M3) supplemented with the equation $\mathit{\text{Det}}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)=0.$ In the case ${\mathit{\lambda}}_{1}=\mathit{\text{Trace}}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)\ne 0$ this system defines coordinates of saddlenode point $\mathit{B}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)$ whose second eigenvalue λ_{2} = 0. The last system supplemented with the equation $\mathit{\text{Trace}}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)=0$ defines an additional degeneration in the system in the point $\mathit{B}\left(\overline{\overline{\mathit{x}}},\overline{\overline{\mathit{z}}}\right)$ where λ_{1} = λ_{2} = 0.
Lemma 1. For a wide range of fixed parameters b, d, s there exist a point $\left(\overline{\mathit{M}},\overline{\mathit{k}}\right)$in the (M, k) parameter plane such that the BogdanovTakens bifurcation of codimension 2 is realized in the model (M1), a = 1 under variations of M and k close to $\left(\overline{\mathit{M}},\overline{\mathit{k}}\right).$The values $\left(\overline{\mathit{M}},\overline{\mathit{k}}\right)$ and coordinates of the equilibrium $\mathit{B}\left(\overline{\overline{\mathit{x}}},\overline{\overline{\mathit{z}}}\right)$where $\phantom{\rule{0.24em}{0ex}}\overline{\overline{\mathit{x}}}=\mathit{x}\left(\overline{\mathit{M}},\overline{\mathit{k}}\right),\overline{\overline{\mathit{z}}}=\mathit{z}\left(\overline{\mathit{M}},\overline{\mathit{k}}\right)$are defined by the system consisting of equations (M3) and equations $\mathit{\text{Det}}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)=0$, $\mathit{Trace}\left(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)\right)=0$(see, e.g. , [41]).
We have found this bifurcation for some reasonable fixed values of the parameters d, b, s with the help of computer package LOCBIF [55]. Using the Lemma and Proposition 1 we prove the following statement.
Proposition 2. (1) System ( M1 ), (3) has equilibria: the saddle O(0,0) and stablenode A(1, 0) for all positive values of the parameters b, d, M, k = 1, s;
(2) there exists positive M* such that

a)
for M < M* the system has only the equilibria O and stable equilibrium A;

b)
for M > M* the system has two more equilibria, a saddle ${\mathit{B}}_{1}\left({\overline{\mathit{x}}}_{1},{\overline{\mathit{z}}}_{1}\right)$ and a stable topological node/spiral
$$\phantom{\rule{0.12em}{0ex}}{\mathit{B}}_{2}\left({\overline{\mathit{x}}}_{2},{\overline{\mathit{z}}}_{2}\right);$$ 
c)
there exist M ^{* *} > M* such that for M > M ^{* *} the spiral $\phantom{\rule{0.12em}{0ex}}{\mathit{B}}_{2}\left({\overline{\mathit{x}}}_{2},{\overline{\mathit{z}}}_{2}\right)$ is placed inside an unstable limit cycle.
The bifurcation diagram of the system under variation of the parameter M is shown in Additional file 1.
Threecomponent model: proof of Statement 1
Let us formulate Statement 1 in more details:

(1)
Equilibrium O(x = y = z = 0) has eigenvalues λ _{1}(O) = 1, λ _{2}(O) =  d, λ _{3}(O) = 1  l; it is unstable for all parameter values;

(2)
If a > 0 then equilibrium A(0,1/a,0) has eigenvalues λ _{1}(A) =  l, λ _{2}(A) =  1, $\phantom{\rule{0.12em}{0ex}}{\mathrm{\lambda}}_{3}\left(\mathrm{A}\right)=\frac{\mathit{da}+\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)}{\mathit{a}};$ it is stable for $\mathit{M}<\frac{\mathit{da}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)}$ and unstable for $\mathit{M}>\frac{\mathit{da}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)}.$
Proof. Jacobian of system (1) is of the form: J(x, y, z) = (a_{i,j})i, j = 1, 2, 3, where
a_{1,1} = 1  l  ax  a(x + y)  bz(1  p(z)), a_{1,2} =  ax + esz, a_{1,3} = esy  bx(1  p(z)) + bxzp_{ z }(z); a_{2,1} = l  ay, a_{2,2} = 1  ay  a(x + y)  b(1  s)z  esz, a_{2,3} =  b(1  s)y  esy,
p_{ z }(z) = 0, if p = const and p_{ z }(z) = k(s  p(z)), if p(z) = (1  s)e^{ kz} + s.
So, the equilibrium O(0, 0, 0) is unstable for both Malthusian and logistic models, the equilibrium $\mathit{A}\left(0,\frac{1}{\mathit{a}},0\right)$ of logistic model is stable for $\mathit{M}<\frac{\mathit{ad}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)}$ and unstable for $\mathit{M}>\frac{\mathit{ad}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)}.$The Statement is proven.
Threecomponent Malthusian model with constant immunity p, p ≥ s > 0
Let us analyze the dynamics of model (1) depending on the parameters l, e. The system has trivial unstable equilibrium O(0,0,0). Let firstly p = s. Then system (1) by changing of variables u = x + y, z = z is reduced to 2component Malthusian system (M1) with respect to u, z  variables. For $\mathit{s}=\mathit{p}<\frac{\mathit{M}}{\mathit{M}+1}$ the orbits {u, z} of this system belongs to the positive quadrant, the nontrivial equilibrium $\overline{\mathit{u}}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)},\overline{\mathit{z}}=\frac{1}{\mathit{b}\left(1\mathit{s}\right)}$ is a center, i.e., it is located inside closed orbits (similar to model (M1)) and trajectories demonstrate periodic oscillations. In variables x(t), y(t), z(t) model (1) also demonstrates periodic oscillations; the model has equilibrium B(x_{ e }, y_{ e }, z_{ e }) where
If p > s then any nontrivial equilibrium (x_{ e }, y_{ e }, z_{ e }) of model (1) is such that the coordinate z_{ e } solves the quadratic equation
and coordinates x_{ e }, y_{ e } can be expressed via z = z_{ e } as follows:
Let $\mathit{l}\left(\mathit{e}\right)=\left(\mathit{M}+\mathit{p}+\mathit{Mp}\right)\left(1+\frac{\left.\mathit{s}(1+\mathit{M}\right)}{\left.\mathit{M}\mathit{s}(1+\mathit{M}\right)}\mathit{e}\right).$ This equation defines a boundary line Γ = (e, l(e)) in the parametric plane (e, l).
Proposition 3.

1)
System (1) with p = const and a = 0 has at most one positive equilibrium B(x _{ e }, y _{ e }, z _{ e }) where x _{ e }, y _{ e }, z _{ e } are defined by formulas (M5), (M6) for s < p and by formula (M4) for s = p;

2)
the system has a single positive equilibrium B (x _{ e }, y _{ e }, z _{ e }) if and only if one of the following conditions holds:

a)
$\mathit{s}<\mathit{p}<\frac{\mathit{M}}{\mathit{M}+1};$

b)
$\mathit{s}<\frac{\mathit{M}}{\mathit{M}+1}<\mathit{p}$ and l > l(e);

3)
in the case a) the equilibrium B is asymptotically stable;

4)
for s = p the system demonstrates periodic oscillations of variables x ( t ), y(t) while z(t) tends to a stable value for a wide domain of initial values (x _{ 0, } y _{0,} z _{0}) close to the equilibrium B.
Proof.
Taking the solutions of quadratic equation (M5) in the form
where D = l^{2}(1  s)^{2} + (p  s + es)^{2}  2l(p(1  s + 2es)  s(1 + e  s + es) we can easily verify that both "branches" z_{1}(l, e), z_{2}(l, e) are real for any positive (e, l) because the expression under the radical is nonnegative. The branches z_{1}(l, e), z_{2}(l, e) are positive both if l < 1 and only z_{1}(l, e) is positive if l > 1.
Analysis of formulas (M5), (M6) shows that only the branch z_{1}(l, e) can define positive coordinates of the equilibrium x_{ e } = x(z_{1}), y_{ e } = y_{ e }(z_{1}). Substituting z_{1}(l, e) into formulas (M6) we obtain that x_{ e }(e, l), y_{ e }(e, l) are positive if the point (e, l) in the parametric plane is placed above the boundary line Γ given by equation
Statements 1 and 2 of the Proposition are proven.
Let us analyze a stability of equilibrium B (x_{ e }, y_{ e }, z_{ e }) of the system. For p = s characteristic polynomial of the system in the point B, whose coordinates are given by (M4), is of the form $\mathit{E}\left({\mathit{\mu}}_{0}\right)=\mathit{Det}(\left(\mathit{J}\left(\mathit{B}\right){\mathit{\mu}}_{0}\mathit{I}\right)=\frac{\left(\mathit{d}+{\mathit{\mu}}_{0}^{2}\right)\left(\mathit{b}\left(\mathit{l}+{\mathit{\mu}}_{0}\right)\left(1\mathit{s}\right)+\mathit{es}\right)}{\mathit{b}\left(1\mathit{s}\right)}$. Thus, two eigenvalues of the point are imaginary, $\phantom{\rule{0.25em}{0ex}}{{\mathit{\mu}}_{0}}^{1,2}=\pm \mathit{i}\sqrt{\mathit{d}},\phantom{\rule{0.25em}{0ex}}$ and the third is negative, $\phantom{\rule{0.25em}{0ex}}{{\mathit{\mu}}_{0}}^{3}=\frac{\mathit{es}+\mathit{bl}\left(1\mathit{s}\right)}{\mathit{b}\left(1\mathit{s}\right)}.$Thus, statement 4 is proven.
We show now that for p > s the point B is a sink, i.e., its eigenvalues have negative real parts (more exactly, one eigenvalue is real negative, and two others are complex with negative real part).
Introduce the parameter α = p  s and write the right hands of (1), a = 0 in the form:
From the condition R(x, y, z) = 0, Q(x, y, z) = 0 we can express x = x_{ e }, y = y_{ e } via z = z_{ e }:
Substituting x = x_{ e }, y = y_{ e } to P(x, y, z) we obtain:
Solving the equation H(z) = 0 we get two roots z_{e1}, z_{e2}; supposing z = z_{0} + αh_{ z } and expanding H(z) in series by α we get the two roots up to o(α) in the form:
Notice, that the root z_{e1} is positive, whereas z_{e2} is positive only for 0 < l < 1. The coordinates x = x_{ e }(α), y = y_{ e }(α) of B are both positive only for z = z_{e1}.
Letting x_{e1} = x_{01} + αh_{1x}, y_{e1} = y_{01} + αh_{1y}, we get
Thus, coordinates of the equilibrium B_{ e } can be presented in the form (x_{ e }, y_{ e }, z_{ e }) = (x_{1e} + αh_{1x}, y_{1e} + αh_{1y}, z_{1e} + αh_{1z}), see formulas (M7), (M8).
We apply the method of small parameter expansion to verify stability of the equilibrium B_{ e }. Characteristic polynomial of the system can be presented in the form
where μ = μ_{0} + αh_{ μ } and
Let now ${\mathit{\mu}}_{0}^{2}=\mathit{d}.$ Substituting this value to Z(μ) and solving equation Z(μ) = 0, we find
and Re(${\mathit{h}}_{\mathit{\mu}})=\frac{\mathit{e}\mathrm{el}\mathit{s}}{2{\left(1\mathit{s}\right)}^{3}}\left(\frac{\mathit{es}\left(\mathit{M}\left(1+\mathit{s}\right)\mathit{s}\right)+\mathit{b}\left(1\mathit{s}\right)\left(\mathit{l}\left(\mathit{M}\left(1+\mathit{s}\right)\mathit{s}\right)\mathit{d}\left(1+\mathit{M}\right)\left(1\mathit{s}\right)\right)}{{\mathit{b}}^{2}\mathit{d}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)}\frac{{\left(1\mathit{s}\right)}^{2}}{\mathit{es}+\mathit{bl}\left(1\mathit{s}\right)}\right)<0$for reasonable parameter values. Thus, equilibrium B_{ e } is asymptotically stable for s < p.
The Proposition is proven.
Threecomponent logistic model with constant immunity p
The logistic version of model (1) with p = s is reduced to 2component logistic system (M1) with respect to variables u = x + y, z. For $\mathit{s}=\mathit{p}<\frac{\mathit{M}}{\mathit{M}+1}$ it has two equilibria, A(0, 1/a) and $\mathit{B}\left(\mathit{u}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)},\mathit{z}=\frac{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\mathit{ad}}{{\mathit{b}}^{2}\left(1\mathit{s}\right)\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)}\right).$ According to Proposition 2, these equilibria are stable in different parameter domains.
In variables x, y, z the equilibrium B(x_{ e }, y_{ e }, z_{ e }) of system (1) has coordinates
Let p be a constant, p > s > 0 . According to Statement 1, the system has trivial equilibrium O(x = 0, y = 0, z = 0), which is unstable for all parameter values, and the equilibrium $\phantom{\rule{0.24em}{0ex}}\mathit{A}\left(\mathit{x}=0,\mathit{y}=\frac{1}{\mathit{a}},\mathit{z}=0\right),$ which is stable if $\mathit{M}<\frac{\mathit{ad}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)}$ and unstable if $\mathit{M}>\frac{\mathit{ad}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)}.$ The system has also a nontrivial equilibrium B whose x, y  coordinates are expressed via z coordinate:
and z coordinate solves the equation:
Denote
then $\overline{\mathit{z}}$ is a root of the equations p = h^{±}(z). Solutions of the latter equation are the points of intersection of the line 0 < p ≤ 1 and the curves h^{±}(z). Two cases with different values of the parameter M are presented in Additional file 2. It demonstrates that at most two positive values of $\overline{\mathit{z}}$ are possible and hence the system may have at most two positive equilibria, whose x, y ‒ coordinates are expressed via $\overline{\mathit{z}}$ by the equations (M9).
Let us define the critical value of the parameter M, ${\mathit{M}}^{\mathit{tc}}=\frac{\mathit{ad}+\mathit{bs}}{\mathit{b}\left(1\mathit{s}\right)}.$
Proposition 4. For $\mathit{s}\le \mathit{p}<\frac{\mathit{M}}{\mathit{M}+1}$and fixed positive values of d, b < 1, l, e system (1) with a = 1 has a single stable equilibrium $\phantom{\rule{0.24em}{0ex}}\mathit{A}\left(\mathit{x}=0,\mathit{y}=\frac{1}{\mathit{a}},\mathit{z}=0\right)$if M < M^{tc}and a single stable equilibrium B(x_{0}, y_{0}, z_{0}) if M > M^{tc}; here the coordinates (x_{0}, y_{0}, z_{0}) are defined by formulas (M10, M11) if s < p and by (M9) if s = p. Transcritical bifurcation "changing of stability" of the points A and B happens as M = M^{tc}.
For p = s system (1) has nontrivial point B, whose coordinates are expressed by (M9); it is easily to verify that ${\mathit{u}}_{\mathit{e}}={\mathit{x}}_{\mathit{e}}+{\mathit{y}}_{\mathit{e}}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)},{\mathit{z}}_{\mathit{e}}=\frac{\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\mathit{ad}}{{\mathit{b}}^{2}\left(1\mathit{s}\right)\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)},$ and 1  a(x_{ e } + y_{ e })  b(1  s)z_{ e } = 0. Hence, z_{ e } > 0 if > M^{tc}; at this condition the equilibrium A looses stability according to Statement 1.
Characteristic equation at the equilibrium B for p = s is of the form:
Accounting the equality 1  a(x_{ e } + y_{ e })  b(1  s)z_{ e } = 0, we obtain from the first term of the last equation: μ_{1} = 1  l  au_{ e }  b(1  s)z_{ e }  esz_{ e } =  l  sz_{ e } < 0, and from the second term ${\mathit{\mu}}_{2,3}=\frac{1}{2}\left(\mathit{au}\pm \sqrt{4\left(\mathit{b}\left(\mathit{M}\left(1\mathit{s}\right)\mathit{s}\right)\mathit{ad}\right){\mathit{u}}_{\mathit{e}}+{\left(\mathit{au}\right)}^{2}}\right).$ For = M^{tc}μ_{2} = 0, μ_{3} =  au < 0. Hence, the transcritical (pitchfork) bifurcation happens with points A, B [41], ch.7. When M > M^{tc} the eigenvalues μ_{2,3} have negative real parts. So, the equilibrium point B is stable for p = s and M > M^{tc}. Due to continuity arguments the point B is also stable for p > s if p is close to s. A general case p > s was verified by computation of eigenvalues of complete characteristic equation in B.
Threecomponent Malthusian model with the immunity p = p(z); proof of Statement 2
The coordinates of nontrivial equilibria are given by formulas (M5) and (M6). Solving equation (M5) with respect to p = p(z), we obtain the equation for zcoordinate:
Evidently, the zcoordinate of possible equilibrium does not depend on the parameter M. Next, h(z) → 1, p(z) → s < 1 as z → ∞, so that two general cases of mutual placing of p(z) and h(z) for l > 1 and l < 1 are possible, which are shown in Additional file 3. Equation (M12) has up to two ("small") positive roots z = z_{ e } if l < 1 (see Additional file 3a) and one ("large") positive root z_{ e } if l > 1 (see Additional file 3b). Equilibrium values x_{ e }(z), y_{ e }(z) defined by the formulas (M6) have different signs for small values of z_{ e } and are both positive for large z_{ e }. Accordingly, the system can have only one positive equilibrium. Eigenvalues of this equilibrium can be computed from the equation Det(J(x_{ e }, y_{ e }, z_{ e })  μI) = 0. Using the LOCBIF software [55], we show that one eigenvalue is real and negative but two other are complex with a positive real part. Thus, this equilibrium is unstable. Statement 2 is proven.
Threecomponent logistic model with the immunity p = p(z)
Proof of Statement 3
Let B_{ M }(x, y, z) be a nontrivial stable equilibrium of model (1), (3) whose coordinates (x(M), y(M), z(M)) depend on the parameter M and satisfy system (2); let P(μ) ≡ Det(J(B)  μI) = 0 be the characteristic polynomial of the system around B_{ M }. The supercritical Hopf bifurcation, corresponding to changing of stability of B_{ M } accompanied by appearance a stable limit cycle happens in the system when for some M a pair of eigenvalues becomes imaginary and certain conditions of nondegeneracy are fulfilled (see, for example, [41]). The Hopf bifurcation is supercritical if the first Lyapunov value becomes negative. The existence of this bifurcation in logistic version of model (1), (3) was verified using LOCBIF [55]. The program numerically finds coordinates of equilibrium with imaginary eigenvalues under variation of parameter M and one more parameter of the model (e.g., e, l, or s) for fixed values of other parameters, checks the sign of the first Lyapunov value and verifies the conditions of nondegeneracy formulated in the Hopf theorem. Applying this software we have found the parameter curves of Hopf bifurcation e_{ H }(M), l_{ H }(M); s_{ H }(M) for fixed values of other parameters of the model (see Additional file 4); it proves the assertions of the Statement.
Reviewers’ reports
Reviewer 1: Sandor Pongor
Comment: The CRISPRCas system is of high theoretical and practical interest. On the practical side it is important for designing genome engineering tools, on the theoretical side, it provides a noteworthy example of Lamarckian evolution. By inserting virusderived spacers into CRISPR repeat cassettes, microbes preserve in their genomes the signature of viruses that attack them, which in turn they pass on to their offspring. The evolution of this system is practically intriguing and has been the subject of several agent based modeling studies. As the authors point out, agent systems can be used to model various levels of complexity, however they do not permit a full and rigorous mathematical analysis of all possible behaviors. Berezovskaya and associates present here a differential equation based model of the evolution of CRISPRCas based immunity. They use LotkaVolterra type models that allow the analysis of Malthusian and logistic regimes and show that the models display complex quasichaotic oscillations. Such complex behaviors have not been found either by experiment or by the previous agent based mutations. The results are well underpinned and clearly discussed. The results are especially interesting example of how unexpected complexity emerges in a seemingly simple system. The text is very well written however the authors may want to check it for minor typos, e.g.:
Page 3 challenge (a virus), "in contrast as opposed" to the random, undirected mutations in the Darwinian evolutionary framework  one of them is superfluous.
Page 4 after eqn 1 host reproduction is Malthusian if a = 0 or logistic if a > 0, and in both cases. – sounds like an unfinished sentence.
Authors’ response: We appreciate these comments. The proposed corrections have been made.
Reviewer 2: Sergei Maslov
Comment: The manuscript addresses a very interesting topic of the emergence of chaotic oscillations in population dynamics of phages and their bacterial hosts. To my knowledge, this topic is relatively unexplored in the literature except for "Bifurcation analysis of bacteria and bacteriophage coexistence in the presence of bacterial debris" by Ira Aviram, Avinoam Rabinovitch [44]. Authors should cite this paper and compare and contrast their results to findings of Aviram et al.
Authors’ response: As far as we can see, the work of Aviram and Rabinovitch is relevant only in a very generic sense. In the revision, we cite this paper along with a more recent publication of the same authors, and comment on the emergence of complex behavior in these models.
Comment: I found the paper very hard to read and understand. In its present form it is unnecessarily heavy on mathematical details and terminology and light on biological insights. I would strongly recommend a *near complete rewriting* of the text of the manuscript that would delegate unnecessary mathematical details and theorems to supplementary materials and explaining the biological consequences of main findings. For example, on page 5 authors refer to their model as "conservative". It is not at all obvious to the majority of even sophisticated computational biology readers that in a conservative system small deviations from the steady state solution do not decay back to the steady state but persist indefinitely. Whenever possible authors should avoid using mathematical jargon and explain in plain English what their parameters/assumptions mean biologically.
Authors’ response: As per this comment and analogous comments of Dr. Kimmel (see below), the entire description of twodimensional (or twocomponent, using the modified terminology), which included most of the mathematical detail, was moved to the Methods. There, we considered it appropriate to formally mathematically define "conservative systems" and other concepts that might be unfamiliar to nonspecialist readers. In the main text, biological interpretations of the parameters and assumptions were provided on several occasions. We believe that these changes made the paper considerably more straightforward, and we appreciate these suggestions of the reviewers. In general, however, one has to face the fact that this is a mathematical biology paper. Further, we beg to disagree that this paper is "light on biological insights". We think that the revealed complex behavior of the coevolving virushost systems is a potentially important biological instant, and the reviewers do not seem to disagree. What is somewhat lacking, are direct connections to experimental results. Such relevant results could come, first, from quantitative measurements of virus diversity and CRISPRCas prevalence in various habitats, and second, from laboratory coevolution experiments. We expect that such results are indeed forthcoming and hope that the present facilitates their interpretation but at this time, there is little data for direct comparison.
Comment: Another example of the same point: on page 5 authors introduce two steady state solutions A and B but do not explain what they mean biologically: phages coexist with bacteria in A but die off in B. And this is just one example of the lack of biological interpretation of mathematical results happening throughout the manuscript.
Authors’ response: This specific statement has been moved to the Methods as per the suggestion of Dr. Maslov and Dr. Kimmel. The interpretation of these steady solutions given by Dr. Maslov is quite correct and thus, we presume, is obvious from the presentation. Especially in the Methods section, further clarifications seem unnecessary. On several other occasions (see below), however, we did include additional biological interpretations.
Comment: On the same page authors cite earlier studies with empirical results confirming that the fraction of immune bacteria p directly depends on the phage population z and not on the bacterial population x. A more detailed discussion of what the empirical data actually say would be beneficial here.
Author’s response: This discussion has been expanded and made more specific.
Comment: In particular, I don’t expect p(T) to instantaneously trace rapidly growing phage population z(T), which seems to be a prerequisite for chaotic behavior reported in this manuscript. What would happen when a delay or time averaging is added to the model? That is to say, what would happen if p(T) is a function of z(Tt_delay) or (in a separate ariant of the model) p(T) is determined by the timeaveraged value of z(T) over T:Tt_average time interval? I would be particularly interested to see if chaotic dynamics would disappear for large enough values of t_delay or t_average. What is a realistic value of t_average compared to a single generation of lytic phage growth?
Authors’ response: These are interesting and important questions that, however, are beyond the scope of the present paper.
Comment: I also request a small yet important change in notation: authors repeatedly refer to threedimensional version of their model which I understood as a model with 3 spatial coordinates. Indeed, special inhomogeneity is known to play an important role in phagebacterial interactions (see e.g. [56, 57]. However, what authors meant is simply a threespecies or threecomponent model with phages, immune hosts, and susceptible hosts. To avoid confusion the term "threedimensional model" while mathematically correct needs to be changed to "three component model" throughout the manuscript.
Authors’ response: We adopted this change in the revision.
Comment: Are figures of plots in Figure 1, 2, 3, 4, 5 necessary? I personally did not learn anything from them. Figure 6, Additional files 1, 2, 3 and 4 nicely illustrate chaotic dynamic of the system. Perhaps they can be collected as multiple panels of just one figure?
Authors’ response: We agree, figures 1, 2, 3, 4, 5have been moved to the additional files.
Additional comments made in the second round of review
Comment: On the other hand, steady state and simple timedependent solutions to LotkaVolterra equations for bacterialphage systems have been studied for quite some time. Some classic references such as [22, 23] have been overlooked and need to be cited in the manuscript.
Authors’ response: We agree, these are relevant references that are cited in the revised version of the present article.
Comment: Population cycles and fixed points in modified LotkaVolterra equations have been also considered in [58, 59]. Would authors predictions of chaotic (or quasichaotic) behavior persist in these systems?
Authors’ response: We do not see how to directly link these models with our approach (at least not without additional, extensive analysis) and therefore currently cannot make such predictions.
Comment: Spatial (see e.g. [56]) and temporal [57] inhomogeneity of the environment is known to play an important role in phagebacterial interactions. How much would it affect chaotic dynamics. These are not idle questions since they go to the heart of the question of how generic is the chaotic behavior reported in the manuscript. If authors believe they are beyond the scope of the current paper, perhaps, these questions should be mentioned in the discussion section as model generalizations/modifications that need to be performed in future studies.
Authors’ response: we fully agree that these are among desirable generalizations of the present approach. It is another matter whether, with the addition of these nonhomogeneities, the model remains tractable. On very general grounds, given that here we have shown that the pseudochaotic oscillations only emerge in a system with certain minimal complexity (distinguishing susceptible and immune hosts is essential), we would expect that such oscillations only become more prominent in even more complex models. However, this is obviously only a conjecture at this point, we cannot be confident before the actual analysis is done.
Comment: Throughout the manuscript authors consider only virulent (lytic) phages. Would there be any interesting modification of predicted dynamical patterns for temperate phages?
Authors’ response: As such, temperate phages, by definition, do not kill the host, and therefore, even if lysogenization is prevented by CRISPRCas, as indeed has been reported [60], this seems to be irrelevant for the modeling approach described here. The situation certainly changes when it comes to prophage induction, against which CRISPRCas protects as well [60]. This case does not appear to be distinguishable from lytic infection within the approximations of the model.
Comment: I appreciated authors following my request and renaming two/three dimensional model into two/three component model. However, on page 11 and several other places in the manuscript old notation is still being used. I recommend authors do global search for "2D", "3D", and " dimensional" in the manuscript.
Authors’ response: This has been taken care of.
Reviewer 3: Marek Kimmel
This paper addresses the issue of coevolution of a virus an and the immune system of a host, taking into account the dynamics of the virus and two types of immune systems: susceptible and resistant. The model is inspired by a type of immune response (CRISPCas) in archaea and some bacteria. The dynamics is summarized by a system of 3 nonlinear ordinary differential equations (ODEs). The system seems to exhibit various dynamical regimes including some that are chaotic. This, according to the authors, provides some analogy to the known examples of the CRISPCas system behavior.
The paper should be reorganized before it is publishable.
Major issues
Comment: 1. The paper is written in a way which makes understanding it very difficult. A large portion of the paper is devoted to models which are inadequate in that they do not include sensitive and resistant immune systems, are therefore limited to two ODEs and cannot exhibit complicated dynamics. To make the paper readable, it should proceed directly to the point. The auxiliary models can be moved to an appendix.
Authors’ response: This reorganization of the manuscript has been implemented as suggested.
Comment: 2. Dynamics of the really interesting 3ODE models is explored mostly numerically, if I understand correctly. In my opinion, more illustrative material might be provided, using the space available after removal of the 2ODE models.
Authors’ response: We carefully considered this suggestion but found that the comparison of Figures 2, 3, 4, 5and 6was highly illustrative of the results for the 3ODE models. The transitions between the outcomes depending on the parameters was made fully explicit in the revised description.
Comment: 3. I am missing a detailed discussion of how the model is similar to experimental observations. What type of data are available? Time series are probably most desirable. Is there any information in the data regarding sensitivity to parameter change and initial conditions, which are characteristic of chaos?
Authors’ response: As pointed out in our response to Dr. Maslov’s comments above, the experimental data are simply not ready for detailed comparison. We would be more than happy to analyze time series or cite an appropriate analysis but such experiments belong in the future.
Comment: 4. Finally, why would this quite generic mathematical description be specific to the CRISPCas systems?
Authors’ response: We never claimed that this description was specific to CRISPRCas. It is inspired by CRISPRCas but is applicable to any system of adaptive immunity with sufficiently long term memory, as pointed out both in the abstract and in the Conclusions.
Detailed remarks
Comment: 1. Page 9. It seems x and y are each composed of a number of different "immunotypes" of host individuals. The structure will be continuous and not twopoint as it is now. What will be the consequences for the dynamics? Will it not be more regular because of smoothing effect of continuity?
Authors’ response: To the best of our understanding, x is just one type. However, y indeed can be represented as numerous "immunotypes" if immunity to different viruses is considered separately. This situation has been explored within the framework of agentbased models [28, 32]. Under the analytic approach used here, continuous distribution of "immunetypes" would inevitably make the model intractable.
Comment: 2. Page 15. Computational analysis usually cannot "show" that an equilibrium is stable. It may be at best consistent with stability.
Authors’ response: The revised version of the manuscript includes a comprehensive test for stability using LOCBIF. Thus, "show" (which does not mean "proven") seems appropriate.
3. "Proposition 4" does not seem to be mathematically demonstrated. So, it is a Conjecture.
Authors response: A sketch of the proof is given in the revision.
References
 1.
Stern A, Sorek R: The phagehost arms race: shaping the evolution of microbes. Bioessays. 2011, 33 (1): 4351.
 2.
Lotka AJ: Elements of Physical Biology. 1925, NY: Williams and Wilkins
 3.
Volterra V: Variazioni e fluttuazioni del numero d’individui in specie animali conviventi. Mem Acad Lincei Roma. 1926, 2: 31113.
 4.
May RM: Stability and Complexity in Model Ecosystems. 2001, Princeton: Princeton University Press
 5.
Volterra V: Le cons sur la Theorie Mathematique de la Lutte popur la Vie. 1931, Paris: Gauthier Villare
 6.
CRISPRCas Systems. RNAmediated Adaptive Immunity in Bacteria and Archaea. Edited by: Barrangou R, Oost J Van der. 2013, Heidelberg: Springer
 7.
Wiedenheft B, Sternberg SH, Doudna JA: RNAguided genetic silencing systems in bacteria and archaea. Nature. 2012, 482 (7385): 331338.
 8.
Marraffini LA, Sontheimer EJ: CRISPR interference: RNAdirected adaptive immunity in bacteria and archaea. Nat Rev Genet. 2010, 11 (3): 181190.
 9.
Makarova KS, Wolf YI, Snir S, Koonin EV: Defense islands in bacterial and archaeal genomes and prediction of novel defense systems. J Bacteriol. 2011, 193 (21): 60396056.
 10.
Koonin EV, Makarova KS: CRISPRCas: Evolution of an RNAbased adaptive immunity system in prokaryotes. RNA Biol. 2013, 10 (5): 679686.
 11.
Sorek R, Lawrence CM, Wiedenheft B: CRISPRmediated adaptive immune systems in bacteria and archaea. Annu Rev Biochem. 2013, 82: 237266.
 12.
Horvath P, Barrangou R: CRISPR/Cas, the immune system of bacteria and archaea. Science. 2010, 327 (5962): 167170.
 13.
Koonin EV, Wolf YI: Is evolution Darwinian or/and Lamarckian?. Biol Direct. 2009, 4: 42
 14.
Koonin EV: The Logic of Chance: The Nature and Origin of Biological Evolution. Upper Saddle River. 2011, Upper Saddle River, NJ: FT press
 15.
Sampson TR, Weiss DS: Exploiting CRISPR/Cas systems for biotechnology. Bioessays. 2014, 36 (1): 3438.
 16.
Ran FA, Hsu PD, Wright J, Agarwala V, Scott DA, Zhang F: Genome engineering using the CRISPRCas9 system. Nat Protoc. 2013, 8 (11): 22812308.
 17.
Gasiunas G, Siksnys V: RNAdependent DNA endonuclease Cas9 of the CRISPR system: Holy Grail of genome editing?. Trends Microbiol. 2013, 21 (11): 562567.
 18.
Mali P, Yang L, Esvelt KM, Aach J, Guell M, Dicarlo JE, Norville JE, Church GM: RNAguided human genome engineering via Cas9. Science. 2013, 339 (6121): 823826.
 19.
Jiang W, Bikard D, Cox D, Zhang F, Marraffini LA: RNAguided editing of bacterial genomes using CRISPRCas systems. Nat Biotechnol. 2013, 31 (3): 233239.
 20.
Cho SW, Kim S, Kim JM, Kim JS: Targeted genome engineering in human cells with the Cas9 RNAguided endonuclease. Nat Biotechnol. 2013, 31 (3): 230232.
 21.
Fitch WM: The variety of human virus evolution. Mol Phylogenet Evol. 1996, 5 (1): 247258.
 22.
Campbell AM: Conditions for the existence of bacteriophage. Evolution. 1960, 15: 153165.
 23.
Levin BR, Stewart FM, Chao L: Resourcelimited growth, competition, and predation: A model and experimental studies with bacteria and bacteriophage. Am Nat. 1977, 111: 324.
 24.
Williams HT: Phageinduced diversification improves host evolvability. BMC Evol Biol. 2013, 13: 17
 25.
Held NL, Herrera A, CadilloQuiroz H, Whitaker RJ: CRISPR associated diversity within a population of Sulfolobus islandicus. PLoS One. 2010, 5 (9): e12988
 26.
Weinberger AD, Sun CL, Pluciński MM, Denef VJ, Thomas BC, Horvath P, Barrangou R, Gilmore MS, Getz WM, Banfield JF: Persisting Lowabundance viral sequences shape microbial CRISPRbased immunity. PLoS Comp Biol. 2012, in press
 27.
Levin BR, Moineau S, Bushman M, Barrangou R: The population and evolutionary dynamics of phage and bacteria with CRISPRmediated immunity. PLoS Genet. 2013, 9 (3): e1003312
 28.
Weinberger AD, Wolf YI, Lobkovsky AE, Gilmore MS, Koonin EV: Viral diversity threshold for adaptive immunity in prokaryotes. MBio. 2012, 3 (6): e00412e00456.
 29.
Childs LM, Held NL, Young MJ, Whitaker RJ, Weitz JS: Multiscale model of CRISPRinduced coevolutionary dynamics: diversification at the interface of Lamarck and Darwin. Evolution. 2012, 66 (7): 20152029.
 30.
Levin BR: Nasty viruses, costly plasmids, population dynamics, and the conditions for establishing and maintaining CRISPRmediated adaptive immunity in bacteria. PLoS Genet. 2010, 6 (10): e1001171
 31.
Gandon S, Vale PF: The evolution of resistance against good and bad infections. J Evol Biol. 2014, 27 (2): 303312.
 32.
Iranzo J, Lobkovsky AE, Wolf YI, Koonin EV: Evolutionary dynamics of the prokaryotic adaptive immunity system CRISPRCas in an explicit ecological context. J Bacteriol. 2013, 195 (17): 38343844.
 33.
Han P, Niestemski LR, Barrick JE, Deem MW: Physical model of the immune response of bacteria against bacteriophage through the adaptive CRISPRCas immune system. Phys Biol. 2013, 10 (2): 025004
 34.
Haerter JO, Sneppen K: Spatial structure and Lamarckian adaptation explain extreme genetic diversity at CRISPR locus. MBio. 2012, 3 (4): e00112e00126.
 35.
Haerter JO, Trusina A, Sneppen K: Targeted bacterial immunity buffers phage diversity. J Virol. 2011, 85 (20): 1055410560.
 36.
Bikard D, Marraffini LA: Innate and adaptive immunity in bacteria: mechanisms of programmed genetic variation to fight bacteriophages. Curr Opin Immunol. 2012, 24 (1): 1520.
 37.
Rechavi O: Guest list or black list: heritable small RNAs as immunogenic memories. Trends Cell Biol. 2013, 24 (4): 212220.
 38.
Rimer J, Cohen IR, Friedman N: Do all creatures possess an acquired immune system of some sort?. Bioessays. 2014, 36 (3): 273281.
 39.
Fineran PC, Charpentier E: Memory of viral infections by CRISPRCas adaptive immune systems: acquisition of new information. Virology. 2012, 434 (2): 202209.
 40.
Westra ER, Swarts DC, Staals RH, Jore MM, Brouns SJ, van der Oost J: The CRISPRs, they are achangin’: how prokaryotes generate adaptive immunity. Annu Rev Genet. 2012, 46: 311339.
 41.
Kuznetsov YA: Elements of applied bifurcation theory. 1995, Vienna: Springer
 42.
Makarova KS, Wolf YI, Koonin EV: Comparative genomics of defense systems in archaea and bacteria. Nucleic Acids Res. 2013, 41 (8): 3604377.
 43.
Makarova KS, Wolf YI, Koonin EV: The basic building blocks and evolution of CRISPRcas systems. Biochem Soc Trans. 2013, 41 (6): 13921400.
 44.
Aviram I, Rabinovitch A: Bifurcation analysis of bacteria and bacteriophage coexistence in the presence of bacterial debris. Commun Nonlinear Sci Numer Simul. 2012, 17: 242254.
 45.
Aviram I, Rabinovitch A: Bacteria and lytic phage coexistence in a chemostat with periodic nutrient supply. Bull Math Biol. 2014, 76 (1): 225244.
 46.
Juliano C, Wang J, Lin H: Uniting germline and stem cells: the function of Piwi proteins and the piRNA pathway in diverse organisms. Annu Rev Genet. 2011, 45: 447469.
 47.
Stuwe E, Toth KF, Aravin AA: Small but sturdy: small RNAs in cellular memory and epigenetics. Genes Dev. 2014, 28 (5): 423431.
 48.
Landry CD, Kandel ER, Rajasethupathy P: New mechanisms in memory storage: piRNAs and epigenetics. Trends Neurosci. 2013, 36 (9): 535542.
 49.
Luteijn MJ, Ketting RF: PIWIinteracting RNAs: from generation to transgenerational epigenetics. Nat Rev Genet. 2013, 14 (8): 523534.
 50.
PaezEspino D, Morovic W, Sun CL, Thomas BC, Ueda K, Stahl B, Barrangou R, Banfield JF: Strong bias in the bacterial CRISPR elements that confer immunity to phage. Nat Commun. 2013, 4: 1430
 51.
Delaney NF, Balenger S, Bonneaud C, Marx CJ, Hill GE, FergusonNoel N, Tsai P, Rodrigo A, Edwards SV: Ultrafast evolution and loss of CRISPRs following a host shift in a novel wildlife pathogen, Mycoplasma gallisepticum. PLoS Genet. 2012, 8 (2): e1002511
 52.
Kuno S, Sako Y, Yoshida T: Diversification of CRISPR within coexisting genotypes in a natural population of the bloomforming cyanobacterium Microcystis aeruginosa. Microbiology. 2014, 160 (Pt 5): 903916.
 53.
Makarova KS, Haft DH, Barrangou R, Brouns SJ, Charpentier E, Horvath P, Moineau S, Mojica FJ, Wolf YI, Yakunin AF, van der Oost J, Koonin EV: Evolution and classification of the CRISPRCas systems. Nat Rev Microbiol. 2011, 9 (6): 467477.
 54.
Bazykin AD: Nonlinear Dynamics of Interacting Populations, Volume 11. 1998, Singapore, New Jersey, London, Hong Kong: World Scientific
 55.
Khibnik AI, Kuznetsov YA, Levitin VV, Nikolaev EV: Continuation techniques and interactive software for bifurcation analysis of ODEs and iterated maps. Physica D. 1993, 62: 360371.
 56.
Heilmann S, Sneppen K, Krishna S: Coexistence of phage and bacteria on the boundary of selforganized refuges. Proc Natl Acad Sci U S A. 2012, 109 (31): 1282812833.
 57.
Maslov S, Sneppen K: Welltemperate phage. arXiv:13081646 [qbioPE]. 2013
 58.
Weitz JS, Dushoff J: Alternative stable states in hostphage dynamics. Theor Ecol. 2008, 1: 1319.
 59.
Wang Z, Goldenfeld N: Fixed points and limit cycles in the population dynamics of lysogenic viruses and their hosts. Phys Rev E Stat Nonlin Soft Matter Phys. 2010, 82 (1 Pt 1): 011918
 60.
Edgar R, Qimron U: The Escherichia coli CRISPR system protects from lambda lysogenization, lysogens, and prophage induction. J Bacteriol. 2010, 192 (23): 62916294.
Acknowledgements
We thank Alex Lobkovsky for helpful discussions. YIW, EVK and GPK are supported through intramural funds of the US Department of Health and Human Services (to the National Library of Medicine).
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
FB and GPK performed the mathematical modeling; YIW, EVK and GPK analyzed the results; GPK and EVK wrote the manuscript that was read and approved by all authors.
Electronic supplementary material
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Cite this article
Berezovskaya, F.S., Wolf, Y.I., Koonin, E.V. et al. Pseudochaotic oscillations in CRISPRvirus coevolution predicted by bifurcation analysis. Biol Direct 9, 13 (2014) doi:10.1186/17456150913
Received
Accepted
Published
DOI
Keywords
 Hopf Bifurcation
 Adaptive Immunity System
 Positive Equilibrium
 Chaotic Oscillation
 Stable Limit Cycle