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

Abstract

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.

Background

The arms races between microbes and viruses preying on them often display rich, complex population dynamics [1]. In principle, the dynamics of virus-microbe interactions is analogous to the classical predator–prey models [2–5] but both microbes and viruses evolve much faster than animals such that virus-host 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 CRISPR-Cas (Clustered Regularly Interspaced Short Palindromic Repeats-CRISPR 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 virus-derived spacers into CRISPR repeat cassettes, thus following the Lamarckian modality of evolution that dramatically accelerates adaptation [13]. The rapid adaptation through the activity of CRISPR-Cas 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 CRISPR-Cas 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 virus-host 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 Lotka-Volterra type models of phage-bacteria 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 long-term persistence of multiple lineages of both the microbial hosts and the viruses [24–27].

Several detailed agent-based models of coevolution between viruses and CRISPR-Cas-carrying hosts have been developed and analyzed [26, 28–35]. The agent-based models incorporate the salient features of the CRISPR-Cas system such as the existence of the CRISPR cassette with virus-derived 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 virus-host systems and predict conditions required for the evolutionary maintenance of the CRISPR-Cas immunity, such as a threshold of viral diversity [28].

Agent-based models provide for the exploration of interactions of arbitrary complexity and naturally incorporate the desired level of granularity (e.g. individual-based or lineage-based models) and the stochasticity of the processes involved. However, such models typically possess a high-dimensional 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 Lotka-Volterra type models of interaction between a host with a heritable adaptive immunity system, such as CRISPR-Cas, 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 CRISPR-Cas system behaviors previously found experimentally and through the agent-based modeling. We explore the full spectrum of possible behaviors of this virus-host system, compare the results with those of a more detailed agent-based model [32], and describe a previously unnoticed regime of quasi-chaotic oscillations.

Results

Three-component CRISPR population dynamics: Malthusian and logistic versions

In order to construct "minimal" analytical models that would qualitatively capture the basic regimes of the CRISPR-Cas system behaviors we analyzed, as a preliminary step, a two-component Volterra-type model which describes the dynamics of virus-host system and consider immunity static within the timescale of the model. Our analysis shows that two-component models display simple dynamical regimes (see Methods for details). However, the CRISPR-Cas 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 CRISPR-Cas 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 non-immune ("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 non-immune, 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 3-component model:

$$\begin{array}{l}\frac{\mathit{\text{dx}}}{\mathit{\text{dt}}}=\mathit{x}(1-\mathit{l}-\mathit{a}\left(\mathit{x}+\mathit{y}\right))-\mathit{\text{bxz}}(1-\mathit{p})+\mathit{\text{esyz}}\equiv \mathit{P}(\mathit{x},\mathit{y},\mathit{z}),\\ \frac{\mathit{\text{dy}}}{\mathit{\text{dt}}}=\mathit{y}+\mathit{lx}-\mathit{\text{ay}}\left(\mathit{x}+\mathit{y}\right)-\mathit{\text{byz}}(1-\mathit{s})-\mathit{\text{esyz}}\equiv \mathit{Q}(\mathit{x},\mathit{y},\mathit{z}),\\ \frac{\mathit{\text{dz}}}{\mathit{\text{dt}}}=\mathit{z}(-\mathit{d}+\mathit{\text{bM}}(\mathit{x}\left(1-\mathit{p}\right)+\mathit{y}(1-\mathit{s}))-\mathit{b}(\mathit{\text{xp}}+\mathit{\text{ys}}))\equiv \mathit{R}(\mathit{x},\mathit{y},\mathit{z}).\end{array}$$

(1)

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:

$$\begin{array}{l}\mathit{P}\left(\mathit{x},\mathit{y},\mathit{z}\right)\equiv \mathit{x}(1-\mathit{l}-\mathit{a}(\mathit{x}+\mathit{y}))-\mathit{\text{bxz}}(1-\mathit{p})+\mathit{\text{esyz}}=0,\\ \mathit{Q}\left(\mathit{x},\mathit{y},\mathit{z}\right)\equiv \mathit{y}+\mathit{\text{lx}}-\mathit{\text{ay}}(\mathit{x}+\mathit{y})-\mathit{\text{byz}}(1-\mathit{s})-\mathit{\text{esyz}}=0,\\ \mathit{R}\left(\mathit{x},\mathit{y},\mathit{z}\right)\equiv \mathit{z}(-\mathit{d}+\mathit{\text{bM}}(\mathit{x}(1-\mathit{p})+\mathit{y}(1-\mathit{s}))-\mathit{b}(\mathit{\text{xp}}+\mathit{\text{ys}}))=0.\end{array}$$

(2)

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 non-immune 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 non-trivial 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 CRISPR-Cas 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 agent-based models of virus-host coevolution suggest that this distinction stems from the fact that hyperthermophiles face lower virus loads and diversity than mesophiles providing for higher efficacy of CRISPR-Cas [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:

$$\mathit{p}\left(\mathit{z}\right)=\left(1-\mathit{s}\right){\mathit{e}}^{-\mathit{\text{kz}}}+\mathit{s}$$

(3)

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 non-trivial equilibrium B _e (x _e , y _e , z _e ) such that x = x _e , y = y _e are expressed via z = z _e :

$${\mathit{x}}_{\mathit{e}}=\frac{-\mathit{d}\left(1-\mathit{b}\left(1-\mathit{s}\right)\mathit{z}\right)}{\mathit{b}\left(\mathit{p}-\mathit{s}\right)\left(1+\mathit{M}-\mathit{\text{bz}}\right)},{\mathit{y}}_{\mathit{e}}=\frac{\mathit{d}\left(1-\mathit{b}\left(1-\mathit{p}\right)\mathit{z}\right)}{\mathit{b}\left(\mathit{p}-\mathit{s}\right)\left(1+\mathit{M}-\mathit{\text{bz}}\right)},$$

(4)

where z = z _e solves the equation

$$\begin{array}{l}\left(1-\mathit{l}\right)+\mathit{b}\left(-2+\mathit{l}+\mathit{p}+\mathit{s}-\mathit{\text{es}}\left(1-\mathit{l}\right)\right)\mathit{z}+{\mathit{b}}^2\left(1-\mathit{p}\right)\\ \times \left(1-\mathit{s}+\mathit{\text{es}}\right){\mathit{z}}^2=0.\end{array}$$

(5)

In particular, we show that for a wide domain of the parameter values, the model demonstrates non-periodic 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 quasi-chaotic 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.

Figure 1

Solutions { x ( t ), y ( t ), z ( t )} and phase portrait of the Malthusian model (1), a = 0 with parameter values l  = 0.9,  e = 0.1. Other parameters are M = 100, d = 1, b = 0.01 , k = 0.5, s = 0.1. Initial values x(0) = 0.07, y(0) = 0.15, z(0) = 22.25 are chosen to be close to the equilibrium.

Let us consider a more realistic version of 3-component 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 non-trivial stable equilibria, stable oscillations, and quasi-chaotic 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 non-trivial equilibrium B(x _e , y _e , z _e ) should satisfy system (2); taking a = 1 for simplicity, we present system (2) in the form:

$$\begin{array}{l}\left(\mathit{x}+\mathit{y}\right)-(\mathit{x}+\mathit{y}){}^2-\mathit{b}(\mathit{x}+\mathit{y})\mathit{z}+\mathit{b}(\mathit{\text{px}}+\mathit{\text{sy}})\mathit{z}=0,\\ -\mathit{d}+\mathit{\text{bM}}\left(\mathit{x}+\mathit{y}\right)-\mathit{b}(\mathit{\text{px}}+\mathit{\text{sy}})(\mathit{M}+1)=0,\\ \mathit{y}+\mathit{\text{lx}}-\mathit{y}\left(\mathit{x}+\mathit{y}\right)-\mathit{b}(1-\mathit{s}+\mathit{\text{es}})\mathit{\text{yz}}=0,\\ \mathit{p}=\mathit{p}\left(\mathit{z}\right)=(1-\mathit{s})\mathit{\text{exp}}(-\mathit{\text{kz}})+\mathit{s}\end{array}$$

(6)

It follows from the first two equations of (6) that

$${\left(\mathit{x}+\mathit{y}\right)}^2-\left(\mathit{x}+\mathit{y}\right)\frac{\left(1+\mathit{M}-\mathit{\text{bz}}\right)}{\mathit{M}+1}+\frac{\mathit{\text{dz}}}{\mathit{M}+1}=0.$$

The solution to the last equation is $\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,

$$\begin{array}{l}0\le {\mathit{x}}_{\mathit{e}}+{\mathit{y}}_{\mathit{e}}<1\text{and}{\mathit{z}}_{\mathit{e}}<\frac{1+\mathit{M}}{\mathit{b}}.\end{array}$$

(7)

It means that all possible non-trivial 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 non-trivial equilibrium B(x _e , y _e , z _e ) such that x = x _e , y =  y _e are expressed via z = z _e

$$\begin{array}{l}\mathit{x}=\frac{-\mathit{d}+\mathit{b}{\mathit{f}}^{\pm }\left(\mathit{z}\right)\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)}{\mathit{b}\left(1+\mathit{M}\right)\left(\mathit{p}-\mathit{s}\right)},\\ \mathit{y}=\frac{\mathit{d}-\mathit{b}{\mathit{f}}^{\pm }\left(\mathit{z}\right)\left(\mathit{M}\left(1-\mathit{p}\right)-\mathit{p}\right)}{\mathit{b}\left(1+\mathit{M}\right)\left(\mathit{p}-\mathit{s}\right)},\text{where}\\ {\mathit{f}}^{\pm }\left(\mathit{z}\right)=\frac{1+\mathit{M}-\mathit{\text{bz}}\pm \sqrt{{\left(1+\mathit{M}-\mathit{\text{bz}}\right)}^2-4\mathit{\text{ad}}\left(1+\mathit{M}\right)\mathit{z}}}{2\mathit{a}\left(1+\mathit{M}\right)}\end{array}$$

(8)

and z -coordinate solves the equation:

$$\begin{array}{l}\left(-1+\mathit{l}+\mathit{a}{\mathit{f}}^{\pm }\right)(\mathit{d}-\mathit{b}{\mathit{f}}^{\pm }(\mathit{M}(1-\mathit{s})-\mathit{s}))\\ +(\mathit{\text{des}}-{\mathit{b}}^2{\mathit{f}}^{\pm }(1-\mathit{p})(\mathit{M}(1-\mathit{s})-\mathit{s})+\mathit{b}(\mathit{d}\left(1-\mathit{p}\right)\\ -\mathit{e}{\mathit{f}}^{\pm }(\mathit{M}(1-\mathit{p})-\mathit{p})\mathit{s}\left)\right)\mathit{z}=0.\end{array}$$

(9)

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^crand 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.

Figure 2

Solutions { x ( t ), y ( t ), z ( t )} and phase portrait of the logistic model (1), (3). Trajectories tend to the stable equilibrium B(0.279, 0.013, 20.24); M = 98.225 < M^cr; k = 0.1, s = 0.2, b = 0.05, l = 0.1, e = 0.5.

Figure 3

 A limit cycle appears in the system (1), (3); M = 98.226 > M^cr.

Figure 4

The limit cycle, phase curve; 5000 < t < 30000, M = 98.226.

Figure 5

Trajectories show quasi-chaotic oscillations, phase curves fill a surface; M = 500.

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 3-dimension system states (see, e.g., [41], ch.5) that there exists such one-to-one 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 non-periodic 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 non-regularity 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, quasi-chaotic oscillations (see Figure 6, right panel).

Figure 6

Left panel: phase curves fill a surface, M = 500; right panel: total host population against the virus population, 1000 < t < 100000.

Discussion

The bifurcation analysis of the virus-host system described here reveals complex, and in particular quasi-chaotic, oscillation regimes that so far have not been previously observed experimentally or in agent-based models of the CRISPR-mediated adaptive immunity [26, 28–35]. The patchy distribution of the CRISPR-Cas systems across the microbial diversity, and even within relatively narrow groups of bacteria [10, 42, 43], implies complex dynamics of virus-host coevolution. Here we show that, in order to detect such non-trivial co-evolutionary 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 two-component models examined here. The key factors for the appearance of quasi-chaotic oscillations are the non-linear dependence of the immunity on the amount of viruses and the three-dimensional 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 bacteriophage-host interaction where the complexity of the system is introduced through inclusion of the time-delayed formation of consumable bacterial debris as a result of cell lysis [44, 45].

Table 1 The analyzed models of the coevolution between viruses and CRISPR-Cas-carrying hosts

It should be emphasized that the three-component 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 CRISPR-Cas 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 quasi-chaotic regime of virus-host coevolution described here might be relevant for these immunity systems as well.

Both the Malthusian and the logistic versions of the three-component model demonstrate quasi-chaotic 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 non-trivial 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 non-regularity of the shapes of the trajectories. When M > > M^cr, the system transits to a regime of quasi-chaotic oscillations.

Thus, the models described here make concrete, quantitative predictions that can be directly tested using an experimental set-up for virus-host coevolution [27, 50].

Conclusions

Comparative genomics as well as direct experiments reveal notable evolutionary instability of the CRISPR-Cas 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 CRISPR-Cas, indicating that this immune system is easily lost and gained via HGT. Rearrangements of the CRISPR-Cas loci are also extremely widespread among microbes [43, 53]. Given this evolutionary plasticity, complex and in particular quasi-chaotic regimes of virus-host 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

Two-component model with p = const

The host-virus dynamics may be considered within the framework of the Volterra-type models

$$\begin{array}{l}\mathit{x}`=\mathit{x}(1-\mathit{\text{ax}}-\mathit{bz}\left(1-\mathit{p}\right))\equiv \mathit{P}(\mathit{x},\mathit{z}),\\ \mathit{z}`=\mathit{z}(-\mathit{d}-\mathit{\text{bxp}}+\mathit{\text{bMx}}\left(1-\mathit{p}\right))\equiv \mathit{Q}(\mathit{x},\mathit{z})\\ \mathit{x}\ge 0,\mathit{z}\ge 0;\mathit{a}\ge 0,\mathit{b},\mathit{d},\mathit{M}>0\end{array}$$

(M1)

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

$$\mathit{B}\left(\overline{\mathit{x}}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}\left(1-\mathit{p}\right)-\mathit{p}\right)},\overline{\mathit{z}}=\frac{1-\mathit{a}\overline{\mathit{x}}}{\mathit{b}\left(1-\mathit{p}\right)}\right).$$

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]).

Two-component 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 non-trivial 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

$$\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)=\left(\begin{array}{cc}0& \mathit{b}\overline{\mathit{x}}(-1+\mathit{p}\left(\overline{\mathit{z}}\right)+\overline{\mathit{z}}{\mathit{p}}_{\mathit{z}}(\overline{\mathit{z}}))\\ -\mathit{b}\overline{\mathit{z}}(-\mathit{M}+\left(1+\mathit{M}\right)\mathit{p}(\overline{\mathit{z}}))& -\mathit{b}\left(1+\mathit{M}\right)\overline{\mathit{x}}\overline{\mathit{z}}{\mathit{p}}_{\mathit{z}}(\overline{\mathit{z}})\end{array}\right),$$

and its determinant and trace are:

$$\begin{array}{l}\mathit{\text{Det}}(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right))={\mathit{b}}^2\overline{\mathit{x}}\overline{\mathit{z}}(\mathit{M}-(1+\mathit{M})\mathit{p}(\overline{\mathit{z}}))(1-\mathit{p}(\overline{\mathit{z}})-\overline{\mathit{z}}{\mathit{p}}_{\mathit{z}}(\overline{\mathit{z}})),\\ \begin{array}{l}\mathit{\text{Trace}}(\mathit{J}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right))=-\mathit{b}\left(1+\mathit{M}\right)\overline{\mathit{x}}\overline{\mathit{z}}{\mathit{p}}_{\mathit{z}}(\overline{\mathit{z}})\end{array}.\end{array}$$

(M2)

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 non-trivial 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}}\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

$$\begin{array}{l}\overline{\mathit{x}}=\frac{\mathit{d}}{\mathit{b}\left(\mathit{M}-\mathit{p}\left(\overline{\mathit{z}}\right)\left(1+\mathit{M}\right)\right)},1-\overline{\mathit{x}}-\mathit{b}\overline{\mathit{z}}\left(1-\mathit{p}\left(\overline{\mathit{z}}\right)\right)=0,\\ \mathit{p}\left(\overline{\mathit{z}}\right)<\mathit{M}/\left(\mathit{M}+1\right).\end{array}$$

(M3)

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:

$$\mathit{J}\left(\mathit{x},\mathit{z}\right)=\left(\begin{array}{cc}1-2\mathit{x}-\mathit{bz}(1-\mathit{p}\left(\mathit{z}\right))& -\mathit{bx}(1-\mathit{p}\left(\mathit{z}\right)-\mathit{z}{\mathit{p}}_{\mathit{z}}(\mathit{z}))\\ -\mathit{bz}(-\mathit{M}+\left(1+\mathit{M}\right)\mathit{p}(\mathit{z}))& -\mathit{d}+\mathit{bMx}(1-\mathit{p}\left(\mathit{z}\right))-\mathit{bxp}(\mathit{z})-\mathit{b}(1+\mathit{M})\mathit{xz}{\mathit{p}}_{\mathit{z}}(\mathit{z})\end{array}\right),$$

$$\begin{array}{l}\mathit{J}\left(0,0\right)=\left(\begin{array}{cc}1& 0\\ 0& -\mathit{d}\end{array}\right),\\ \mathit{J}\left(1,0\right)=\left(\begin{array}{cc}-1& -\mathit{b}(1-\mathit{p}\left(0\right))\\ 0& -\mathit{d}+\mathit{\text{bM}}(1-\mathit{p}\left(0\right))-\mathit{\text{bp}}(0)\end{array}\right)=\left(\begin{array}{cc}-1& 0\\ 0& -\mathit{d}-\mathit{b}\end{array}\right).\end{array}$$

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 saddle-node point $\mathit{B}\left(\overline{\mathit{x}},\overline{\mathit{z}}\right)$ whose second eigenvalue λ₂ = 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 λ₁ = λ₂ = 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 Bogdanov-Takens bifurcation of co-dimension 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 $\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 stable-node 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

  $${\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 ${\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.

Three-component model: proof of Statement 1

Let us formulate Statement 1 in more details:

• (1) Equilibrium O(x = y = z = 0) has eigenvalues λ ₁(O) = 1, λ ₂(O) = - d, λ ₃(O) = 1 - l; it is unstable for all parameter values;

• (2) If a > 0 then equilibrium A(0,1/a,0) has eigenvalues λ ₁(A) = - l, λ ₂(A) = - 1, ${\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,

$${\mathit{a}}_{3,1}=\mathit{bMz}\left(1-\mathit{p}\left(\mathit{z}\right)\right)-\mathit{bzp}\left(\mathit{z}\right),{\mathit{a}}_{3,2}=\mathit{bM}\left(1-\mathit{s}\right)\mathit{z}-\mathit{bsz},$$

$${\mathit{a}}_{3,3}=-\mathit{d}+\mathit{b}\left(\mathit{M}\left(\mathit{x}+\mathit{y}\right)-\left(\mathit{M}+1\right)\left(\mathit{p}\left(\mathit{z}\right)\mathit{x}+\mathit{sy}\right)\right)-\mathit{bxz}{\mathit{p}}_{\mathit{z}}\left(\mathit{z}\right)\left(1+\mathit{M}\right).$$

p _z (z) = 0, if p = const and p _z (z) = k(s - p(z)), if p(z) = (1 - s)e^- kz + s.

$$\begin{array}{l}\mathit{J}\left(0,0,0\right)=\left(\begin{array}{ccc}1-\mathit{l}& 0& 0\\ \mathit{l}& 1& 0\\ 0& 0& -\mathit{d}\end{array}\right),\\ \mathit{J}\left(0,1/\mathit{a},0\right)=\left(\begin{array}{ccc}-\mathit{l}& 0& \mathit{bes}/\mathit{a}\\ -1+\mathit{l}& -1& \mathit{b}\left(-1+\mathit{s}-\mathit{es}\right)/\mathit{a}\\ 0& 0& -(\mathit{ad}-\mathit{b}\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)/\mathit{a}\end{array}\right).\end{array}$$

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.

Three-component 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 2-component 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 non-trivial 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

$$\begin{array}{l}{\mathit{x}}_{\mathit{e}}=\frac{\mathit{des}}{\mathit{b}\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)\left(\mathit{bl}\left(1-\mathit{s}\right)+\mathit{es}\right)},\\ {\mathit{y}}_{\mathit{e}}=\frac{\mathit{dl}\left(1-\mathit{s}\right)}{\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)\left(\mathit{bl}\left(1-\mathit{s}\right)+\mathit{es}\right)},\\ {\mathit{z}}_{\mathit{e}}=\frac{1}{\mathit{b}\left(1-\mathit{s}\right)}\end{array}$$

(M4)

If p > s then any non-trivial equilibrium (x _e , y _e , z _e ) of model (1) is such that the coordinate z _e solves the quadratic equation

$$\left(1-\mathit{l}\right)+\mathit{b}\left(-2+\mathit{l}+\mathit{p}+\mathit{s}-\mathit{es}\left(1-\mathit{l}\right)\right)\mathit{z}+{\mathit{b}}^2\left(1-\mathit{p}\right)\left(1-\mathit{s}+\mathit{es}\right){\mathit{z}}^2=0,$$

(M5)

and coordinates x _e , y _e can be expressed via z = z _e as follows:

$${\mathit{x}}_{\mathit{e}}=\frac{-\mathit{d}\left(1-\mathit{b}\left(1-\mathit{s}\right)\mathit{z}\right)}{\mathit{b}\left(\mathit{p}-\mathit{s}\right)\left(1+\mathit{M}-\mathit{bz}\right)},{\mathit{y}}_{\mathit{e}}=\frac{\mathit{d}\left(1-\mathit{b}\left(1-\mathit{p}\right)\mathit{z}\right)}{\mathit{b}\left(\mathit{p}-\mathit{s}\right)\left(1+\mathit{M}-\mathit{bz}\right)}.$$

(M6)

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 ₀) close to the equilibrium B.

Proof.

Taking the solutions of quadratic equation (M5) in the form

$$\begin{array}{l}\mathit{z}={\mathit{z}}_1\left(\mathit{l},\mathit{e}\right)=\frac{2-\mathit{l}-\left(\mathit{p}+\mathit{s}\right)+\mathit{es}(1-\mathit{l})+\sqrt{\mathit{D}}}{2\mathit{b}\left(1-\mathit{p}\right)\left(1-\mathit{s}+\mathit{es}\right)},\\ \mathit{z}={\mathit{z}}_2\left(\mathit{l},\mathit{e}\right)=\frac{2-\mathit{l}-\left(\mathit{p}+\mathit{s}\right)+\mathit{es}(1-\mathit{l})+\sqrt{\mathit{D}}}{2\mathit{b}\left(1-\mathit{p}\right)\left(1-\mathit{s}+\mathit{es}\right)}\end{array}$$

where D = l²(1 - s)² + (p - s + es)² - 2l(p(1 - s + 2es) - s(1 + e - s + es) we can easily verify that both "branches" z₁(l, e), z₂(l, e) are real for any positive (e, l) because the expression under the radical is non-negative. The branches z₁(l, e), z₂(l, e) are positive both if  l < 1 and only z₁(l, e) is positive if  l > 1.

Analysis of formulas (M5), (M6) shows that only the branch z₁(l, e) can define positive coordinates of the equilibrium x _e  = x(z₁), y _e  = y _e (z₁). Substituting z₁(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

$$\begin{array}{l}\left(-1+\mathit{p}\right)\left(1+\left(-1+\mathit{e}\right)\mathit{s}\right)(\mathit{l}\left(\mathit{M}\left(-1+\mathit{s}\right)+\mathit{s}\right)\\ +\left(\mathit{M}\left(-1+\mathit{p}\right)+\mathit{p}\right)\left(\left(-1+\mathit{e}\right)\mathit{s}+\mathit{M}\left(1+\left(-1+\mathit{e}\right)\mathit{s}\right)\right)=0.\end{array}$$

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, ${{\mathit{\mu }}_0}^{1,2}=\pm \mathit{i}\sqrt{\mathit{d}},$ and the third is negative, ${{\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:

$$\begin{array}{l}\mathit{P}\left(\mathit{x},\mathit{y},\mathit{z}\right)=\mathit{x}-\mathit{lx}-\mathit{b}(1-\mathit{\alpha }-\mathit{s})\mathit{xz}+\mathit{esyz},\\ \mathit{Q}\left(\mathit{x},\mathit{y},\mathit{z}\right)=\mathit{lx}+\mathit{y}-\mathit{b}(1-\mathit{s})\mathit{yz}-\mathit{esyz},\\ \mathit{R}\left(\mathit{x},\mathit{y},\mathit{z}\right)=\mathit{z}(-\mathit{d}+\mathit{bM}((1-\mathit{\alpha }-\mathit{s})\mathit{x}+(1-\mathit{s})\mathit{y})\\ -\mathit{b}((\mathit{\alpha }+\mathit{s})\mathit{x}+\mathit{sy})).\end{array}$$

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 :

$$\begin{array}{l}{\mathit{x}}_{\mathit{e}}=\frac{-(\mathit{d}\left(1-\mathit{b}\left(1-\mathit{s}\right)\mathit{z}+\mathit{esz}\right)}{\mathit{b}(\mathit{\alpha }\left(1+\mathit{M}\right)\left(1-\mathit{b}\left(1-\mathit{s}\right)\mathit{z}+\mathit{esz}\right)+\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)\left(-1+\mathit{l}+\mathit{esz}+\mathit{bz}\left(1-\mathit{s}\right)\right)}.\\ {\mathit{y}}_{\mathit{e}}=\frac{\mathit{dl}}{\mathit{b}(\mathit{\alpha }\left(1+\mathit{M}\right)\left(1-\mathit{b}\left(1-\mathit{s}\right)\mathit{z}+\mathit{esz}\right)+\left(\mathit{M}\left(1-\mathit{s}\right)+\mathit{s}\right)\left(-1+\mathit{l}+\mathit{esz}+\mathit{bz}\left(1-\mathit{s}\right)\right)}\end{array}$$

Substituting x = x _e , y = y _e to P(x, y, z) we obtain:

$$\mathit{P}\left({\mathit{x}}_{\mathit{e}},{\mathit{y}}_{\mathit{e}},\mathit{z}\right)\equiv \mathit{H}\left(\mathit{z}\right)=\frac{\mathit{d}\left(\mathit{l}\left(1-\mathit{b}\left(1-\mathit{s}\right)\mathit{z}\right)-\left(1-\mathit{b}\left(1-\mathit{s}\right)\mathit{z}-\mathit{esz}\right)\left(1-\mathit{b}\left(1-\mathit{s}-\mathit{\alpha }\right)\mathit{z}\right)\right)}{\mathit{b}(\mathit{\alpha }\left(1+\mathit{M}\right)\left(1-\mathit{b}\left(1-\mathit{s}\right)\mathit{z}+\mathit{esz}\right)+\left(\mathit{M}\left(1-\mathit{s}\right)+\mathit{s}\right)\left(-1+\mathit{l}+\mathit{esz}+\mathit{bz}\left(1-\mathit{s}\right)\right)}.$$

Solving the equation H(z) = 0 we get two roots z_e1, z_e2; supposing z = z₀ + αh _z and expanding H(z) in series by α we get the two roots up to o(α)  in the form:

$$\begin{array}{ll}{\mathit{z}}_{\mathit{e}1}={\mathit{z}}_{01}+\mathit{\alpha }{\mathit{h}}_{1\mathit{z}},& {\mathit{z}}_{\mathit{e}2}={\mathit{z}}_{02}+\mathit{\alpha }{\mathit{h}}_{2\mathit{z}},\\ {\mathit{z}}_{01}=\frac{1}{\mathit{b}\left(1-\mathit{s}\right)},& {\mathit{h}}_{1\mathit{z}}=\frac{\mathit{es}}{\mathit{b}{\left(1-\mathit{s}\right)}^2\left(\mathit{bl}\left(1-\mathit{s}\right)+\mathit{es}\right)};\\ {\mathit{z}}_{02}=\frac{1-\mathit{l}}{\mathit{b}\left(1-\mathit{s}\right)+\mathit{es}},& {\mathit{h}}_{2\mathit{z}}=\frac{\mathit{bl}\left(1-\mathit{l}\right)}{\left(\mathit{b}\left(1-\mathit{s}\right)+\mathit{es}\left)\right(\mathit{bl}\left(1-\mathit{s}\right)+\mathit{es}\right)}.\end{array}$$

(M7)

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₀₁ + αh_1x, y_e1 = y₀₁ + αh_1y, we get

$$\begin{array}{l}{\mathit{x}}_{01}=\frac{\mathit{des}}{\mathit{b}\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)\left(\mathit{bl}\left(1-\mathit{s}\right)+\mathit{es}\right)},\\ {\mathit{h}}_{1\mathit{x}}=\frac{-\mathit{des}\left({\mathit{e}}^2{\mathit{s}}^2\left(1+\mathit{M}\right)+{\mathit{b}}^2\mathit{el}\left(1-\mathit{s}\right)\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)+\mathit{bel}\left(1+2\mathit{M}\left(1-\mathit{s}\right)\mathit{s}-2{\mathit{s}}^2\right)\right)}{\mathit{b}{\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)}^2{\left(\mathit{bl}\left(1-\mathit{s}\right)+\mathit{es}\right)}^3};\\ {\mathit{y}}_{01}=\frac{\mathit{dl}\left(1-\mathit{s}\right)}{\mathit{b}\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)\left(\mathit{bl}\left(1-\mathit{s}\right)+\mathit{es}\right))},\\ {\mathit{h}}_{1\mathit{y}}=\frac{\mathit{dels}\left(\mathit{es}-\mathit{b}\left(1-\mathit{s}\right)\left(\mathit{M}-\mathit{el}\left(1+\mathit{M}\right)\left(1-\mathit{s}\right)+\mathit{s}+\mathit{Ms}\right)\right)}{\mathit{b}{\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)}^2{\left(\mathit{bl}\left(1-\mathit{s}\right)+\mathit{es}\right)}^3}.\end{array}$$

(M8)

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

$$\begin{array}{ll}\mathit{E}\left(\mathit{\mu }\right)=\mathit{Det}(\left(\mathit{J}\left({\mathit{B}}_{\mathit{e}}\right)-\mathit{\mu 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)}\\ -\frac{\mathit{\alpha }}{\mathit{b}{\left(1-\mathit{s}\right)}^2\left(\mathit{M}\left(1-\mathit{s}\right)+\mathit{s}\right)\mathit{bl}\left(1-\mathit{s}\right)+\mathit{es})}\mathit{Z}\left(\mathit{\mu }\right)\end{array}$$

where μ = μ₀ + αh _μ and

$$\begin{array}{ll}\mathit{z}\left(\mathit{\mu }\right)& =(-({\mathit{b}}^2\mathit{l}(\mathit{d}(1+{\mathit{h}}_{\mathit{\mu }}(1-\mathit{s}))+{\mathit{\mu }}_0\left({\mathit{\mu }}_0\left(1-3{\mathit{h}}_{\mathit{\mu }}\left(1-\mathit{s}\right)\right)\right.\\ \left.+2\mathit{l}{\mathit{h}}_{\mathit{\mu }}\left(1-\mathit{s}\right)\right)(((1-\mathit{s}){}^2\mathit{M}(1-\mathit{s})-\mathit{s})\\ -\mathit{e}{\mathit{s}}^2(\mathit{d}+{\mathit{\mu }}_0({\mathit{\mu }}_0+2{\mathit{h}}_{\mathit{\mu }}\left(1-\mathit{s}\right)))((1-\mathit{s}){}^2(\mathit{M}(1-\mathit{s})-\mathit{s})\\ +\mathit{be}(1-\mathit{s})\mathit{s}\left({\mathit{\mu }}_0\right({\mathit{\mu }}_0(-1-3{\mathit{h}}_{\mathit{\mu }}\left(1-\mathit{s}\right))\\ -4\mathit{l}{\mathit{h}}_{\mathit{\mu }}(1-\mathit{s}))(\mathit{M}(-1+\mathit{s})+\mathit{s})\\ +\mathit{d}\left(\left(-1-{\mathit{h}}_{\mathit{\mu }}\left(1-\mathit{s}\right)\right)\right(\mathit{M}(-1+\mathit{s})+\mathit{s})+\mathit{l}(\mathit{M}\left(-1+\mathit{s}\right)\\ -{\mathit{\mu }}_0(1+\mathit{M})(-1+\mathit{s})+\mathit{s}\left)\right)\left)\right)).\end{array}$$

Let now ${\mathit{\mu }}_0^2=-\mathit{d}.$ Substituting this value to Z(μ) and solving equation Z(μ) = 0, we find

$${\mathit{h}}_{\mathit{\mu }}=\frac{\left(\mathit{bdels}\left(\mathit{M}+{\mathit{\mu }}_0\left(1+\mathit{M}\right)\left(-1+\mathit{s}\right)-\mathit{s}-\mathit{Ms}\right)\right)}{\left(2\left(\mathit{M}\left(-1+\mathit{s}\right)+\mathit{s}\right)\left(\mathit{bl}\left(-1+\mathit{s}\right)-\mathit{es}\right)\left(\mathit{b}\left(\mathit{d}-\mathit{l}{\mathit{\mu }}_0\right)\left(-1+\mathit{s}\right)+\mathit{es}{\mathit{\mu }}_0\right)\right)}$$

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.

Three-component logistic model with constant immunity p

The logistic version of model (1) with p = s is reduced to 2-component 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

$$\begin{array}{l}{\mathit{x}}_{\mathit{e}}=\frac{\mathit{des}\left(\mathit{b}\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)-\mathit{ad}\right)}{\mathit{b}\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)\left({\mathit{b}}^2\mathit{l}\left(1-\mathit{s}\right)\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)+\mathit{es}\left(\mathit{b}\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)-\mathit{ad}\right)\right)},\\ {\mathit{y}}_{\mathit{e}}=\frac{\mathit{bdl}\left(1-\mathit{s}\right)}{\mathit{b}\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)\left({\mathit{b}}^2\mathit{l}\left(1-\mathit{s}\right)\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)+\mathit{es}\left(\mathit{b}\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)-\mathit{ad}\right)\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)}.\end{array}$$

(M9)

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 $\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 non-trivial equilibrium B whose x, y - coordinates are expressed via z -coordinate:

$$\begin{array}{l}\mathit{x}=\frac{-\mathit{d}+\mathit{b}{\mathit{f}}^{\pm }\left(\mathit{z}\right)\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)}{\mathit{b}\left(1+\mathit{M}\right)\left(\mathit{p}-\mathit{s}\right)},\\ \mathit{y}=\frac{\mathit{d}-\mathit{b}{\mathit{f}}^{\pm }\left(\mathit{z}\right)\left(\mathit{M}\left(1-\mathit{p}\right)-\mathit{p}\right)}{\mathit{b}\left(1+\mathit{M}\right)\left(\mathit{p}-\mathit{s}\right)},\text{where}\\ {\mathit{f}}^{\pm }\left(\mathit{z}\right)=\frac{1+\mathit{M}-\mathit{bz}\pm \sqrt{{\left(1+\mathit{M}-\mathit{bz}\right)}^2-4\mathit{ad}\left(1+\mathit{M}\right)\mathit{z}}}{2\mathit{a}\left(1+\mathit{M}\right)}\end{array}$$

(M10)

and z -coordinate solves the equation:

$$\begin{array}{l}\left(-1+\mathit{l}+\mathit{a}{\mathit{f}}^{\pm }\right)(\mathit{d}-\mathit{b}{\mathit{f}}^{\pm }(\mathit{M}(1-\mathit{s})-\mathit{s}))+(\mathit{des}-{\mathit{b}}^2{\mathit{f}}^{\pm }(1-\mathit{p})(\mathit{M}(1-\mathit{s})-\mathit{s})\\ +\mathit{b}(\mathit{d}\left(1-\mathit{p}\right)-\mathit{e}{\mathit{f}}^{\pm }(\mathit{M}(1-\mathit{p})-\mathit{p})\mathit{s}))\mathit{z}=0.\end{array}$$

(M11)

Denote

$$\begin{array}{l}{\mathit{h}}^{\pm }\left(\mathit{z}\right)\equiv \frac{\left(\mathit{d}\left(\mathit{l}-1+\mathit{a}{\mathit{f}}^{\pm }\right)+\mathit{b}{\mathit{f}}^{\pm }\left(\mathit{M}\left(1-\mathit{a}{\mathit{f}}^{\pm }-\mathit{l}\left(1-\mathit{s}\right)\right)+\mathit{ls}\right)\right)+\left(\mathit{d}-\mathit{b}{\mathit{f}}^{\pm }\mathit{M}\right)\left(\mathit{b}\left(1-\mathit{s}\right)+\mathit{es}\right)\mathit{z}}{\mathit{b}{\mathit{f}}^{\pm }\left(1+\mathit{M}\right)\left(1-\mathit{a}{\mathit{f}}^{\pm }-\left(\mathit{b}\left(1-\mathit{s}\right)+\mathit{es}\right)\mathit{z}\right)}\end{array}$$

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 $\mathit{A}\left(\mathit{x}=0,\mathit{y}=\frac{1}{\mathit{a}},\mathit{z}=0\right)$if M < M^tcand a single stable equilibrium B(x₀, y₀, z₀) if M > M^tc; here the coordinates (x₀, y₀, z₀) 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:

$$\mathit{Det}(\left(\mathit{J}\left(\mathit{B}\right)-\mathit{\mu I}\right)=-\left(-1+\mathit{l}+\mathit{\mu }+\mathit{a}\left({\mathit{x}}_{\mathit{e}}+{\mathit{y}}_e\right)+\mathit{b}\left(1-\mathit{s}\right){\mathit{z}}_e+{\mathit{sz}}_e\right)*$$

$$\left({\mathit{\mu }}^2+{\mathit{b}}^2\left(1-\mathit{s}\right)\left(\mathit{M}\left(1-\mathit{s}\right)-\mathit{s}\right)\left({\mathit{x}}_{\mathit{e}}+{\mathit{y}}_{\mathit{e}}\right){\mathit{z}}_{\mathit{e}}+\mathit{\mu }\left(-1+2\mathit{a}\left({\mathit{x}}_{\mathit{e}}+{\mathit{y}}_{\mathit{e}}\right)+{\mathit{bz}}_{\mathit{e}}\left(1-\mathit{s}\right)\right)\right)=0.$$

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 - 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μ₂ = 0, μ₃ = - 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.

Three-component Malthusian model with the immunity p = p(z); proof of Statement 2

The coordinates of non-trivial equilibria are given by formulas (M5) and (M6). Solving equation (M5) with respect to p = p(z), we obtain the equation for z-coordinate:

$$\begin{array}{ll}\mathit{p}\left(\mathit{z}\right)\equiv \left(1-\mathit{s}\right){\mathit{e}}^{-\mathit{kz}}+\mathit{s}& =\frac{-\left(1-\mathit{l}\right)+\left(\mathit{b}\left(2-\mathit{l}\left(1-\mathit{s}\right)-\mathit{s}\right)+\mathit{es}\right)\mathit{z}-\mathit{b}\left(\mathit{b}\left(1-\mathit{s}\right)+\mathit{es}\right){\mathit{z}}^2}{\mathit{bz}\left(1-\left(\mathit{b}\left(1-\mathit{s}\right)+\mathit{es}\right)\mathit{z}\right)}\\ \equiv \mathit{h}\left(\mathit{z}\right)\end{array}$$

(M12)

Evidently, the z-coordinate 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.

Three-component logistic model with the immunity p = p(z)

Proof of Statement 3

Let B _M (x, y, z)  be a non-trivial 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 non-degeneracy 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 non-degeneracy 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 CRISPR-Cas 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 virus-derived 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 CRISPR-Cas based immunity. They use Lotka-Volterra type models that allow the analysis of Malthusian and logistic regimes and show that the models display complex quasi-chaotic 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 two-dimensional (or two-component, 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 non-specialist 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 co-evolving virus-host 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 CRISPR-Cas prevalence in various habitats, and second, from laboratory co-evolution 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 co-exist 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(T-t_delay) or (in a separate ariant of the model) p(T) is determined by the time-averaged value of z(T) over T:T-t_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 three-dimensional 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 phage-bacterial interactions (see e.g. [56, 57]. However, what authors meant is simply a three-species or three-component model with phages, immune hosts, and susceptible hosts. To avoid confusion the term "three-dimensional 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 time-dependent solutions to Lotka-Volterra equations for bacterial-phage 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 Lotka-Volterra equations have been also considered in [58, 59]. Would authors predictions of chaotic (or quasi-chaotic) 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 phage-bacterial 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 non-homogeneities, the model remains tractable. On very general grounds, given that here we have shown that the pseudo-chaotic 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 CRISPR-Cas, 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 CRISPR-Cas 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 co-evolution 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 (CRISP-Cas) 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 CRISP-Cas 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 3-ODE 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 2-ODE 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 3-ODE 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 CRISP-Cas systems?

Authors’ response: We never claimed that this description was specific to CRISPR-Cas. It is inspired by CRISPR-Cas 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 "immuno-types" of host individuals. The structure will be continuous and not two-point 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 "immuno-types" if immunity to different viruses is considered separately. This situation has been explored within the framework of agent-based models [28, 32]. Under the analytic approach used here, continuous distribution of "immune-types" 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.