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

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.

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.

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.

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

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.

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.

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.

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.

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

Two-component model with p = const

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.

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

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;

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

Proof. Jacobian of system (1) is of the form: J(x, y, z) = (a_i,j)i, j = 1, 2, 3, where

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

Proof.

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.

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

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.

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

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.

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 ^tc and 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.

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

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.

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.