- Open Access
Pseudo-chaotic oscillations in CRISPR-virus coevolution predicted by bifurcation analysis
© Berezovskaya et al.; licensee BioMed Central Ltd. 2014
- Received: 2 April 2014
- Accepted: 26 June 2014
- Published: 2 July 2014
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.
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.
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.
This manuscript was reviewed by Sandor Pongor, Sergei Maslov and Marek Kimmel. For the complete reports, go to the Reviewers’ Reports section.
- Hopf Bifurcation
- Adaptive Immunity System
- Positive Equilibrium
- Chaotic Oscillation
- Stable Limit Cycle
The arms races between microbes and viruses preying on them often display rich, complex population dynamics . 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 . 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 . 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 . 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 . 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 .
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 , and describe a previously unnoticed regime of quasi-chaotic oscillations.
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 . 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.
Equilibrium O(x = y = z = 0) is unstable for all parameter values;
If a > 0 then equilibrium A(0,1/a,0) is stable for and unstable for
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 .
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 the coordinates (x e , y e , z e ) of this equilibrium satisfy (4), (5).
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.
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., , 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.
The analyzed models of the coevolution between viruses and CRISPR-Cas-carrying hosts
p = const
Center M > p/(1-p)
Damped oscillations M > p/(1-p)
p = p ( z )
No stable equilibrium z > 0
Damped oscillations M > M* (finite or infinite basin of attraction)
p = const
Damped oscillations M > p/(1-p) or s/(1-s) < M < p/(1-p) and l > l(e)
Stable equilibrium for
p = p ( z )
Quasi-chaotic oscillations at all M
Damped oscillations at M < M cr ; periodic oscillations at M ≥ M cr ; quasi-chaotic oscillations at M > > M cr
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 . 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 . 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.
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.
Two-component model with p = const
If a > 0 (logistic case), then for constant system (M1) has a saddle in the origin, equilibria A(1/a, 0) and 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., ).
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 has to be a root of the equation 1 - bz(1 - p(z)) = 0 and The point is positive only if
Proposition 1. If p z (z) < 0, then the Malthusian model has only one non-trivial equilibrium 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 can intersect only once. Thus, equation 1 - bz(1 - p(z)) = 0 has only one root , and the Malthusian model has only one non-trivial equilibrium Next, for positive if p z (z) < 0. Thus 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 In the case this system defines coordinates of saddle-node point whose second eigenvalue λ2 = 0. The last system supplemented with the equation defines an additional degeneration in the system in the point where λ1 = λ2 = 0.
Lemma 1. For a wide range of fixed parameters b, d, s there exist a point 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 The values and coordinates of the equilibrium where are defined by the system consisting of equations (M3) and equations , (see, e.g. , ).
We have found this bifurcation for some reasonable fixed values of the parameters d, b, s with the help of computer package LOCBIF . 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;
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 and a stable topological node/spiral
there exist M * * > M* such that for M > M * * the spiral 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
Equilibrium O(x = y = z = 0) has eigenvalues λ 1(O) = 1, λ 2(O) = - d, λ 3(O) = 1 - l; it is unstable for all parameter values;
If a > 0 then equilibrium A(0,1/a,0) has eigenvalues λ 1(A) = - l, λ 2(A) = - 1, it is stable for and unstable for
Proof. Jacobian of system (1) is of the form: J(x, y, z) = (ai,j)i, j = 1, 2, 3, where
So, the equilibrium O(0, 0, 0) is unstable for both Malthusian and logistic models, the equilibrium of logistic model is stable for and unstable for The Statement is proven.
Three-component Malthusian model with constant immunity p, p ≥ s > 0
Let This equation defines a boundary line Γ = (e, l(e)) in the parametric plane (e, l).
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;
the system has a single positive equilibrium B (x e , y e , z e ) if and only if one of the following conditions holds:
and l > l(e);
in the case a) the equilibrium B is asymptotically stable;
for s = p the system demonstrates periodic oscillations of variables x ( t ), y(t) while z(t) tends to a stable value for a wide domain of initial values (x 0, y 0, z 0) close to the equilibrium B.
where D = l2(1 - s)2 + (p - s + es)2 - 2l(p(1 - s + 2es) - s(1 + e - s + es) we can easily verify that both "branches" z1(l, e), z2(l, e) are real for any positive (e, l) because the expression under the radical is non-negative. The branches z1(l, e), z2(l, e) are positive both if l < 1 and only z1(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 . Thus, two eigenvalues of the point are imaginary, and the third is negative, 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 ze1 is positive, whereas ze2 is positive only for 0 < l < 1. The coordinates x = x e (α), y = y e (α) of B are both positive only for z = ze1.
Thus, coordinates of the equilibrium B e can be presented in the form (x e , y e , z e ) = (x1e + αh1x, y1e + αh1y, z1e + αh1z), see formulas (M7), (M8).
and Re(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 it has two equilibria, A(0, 1/a) and According to Proposition 2, these equilibria are stable in different parameter domains.
then 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 are possible and hence the system may have at most two positive equilibria, whose x, y ‒ coordinates are expressed via by the equations (M9).
Let us define the critical value of the parameter M,
Proposition 4. For and fixed positive values of d, b < 1, l, e system (1) with a = 1 has a single stable equilibrium if M < M tc and a single stable equilibrium B(x0, y0, z0) if M > M tc ; here the coordinates (x0, y0, z0) 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 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 = 1 - l - au e - b(1 - s)z e - esz e = - l - sz e < 0, and from the second term For = M tc μ2 = 0, μ3 = - au < 0. Hence, the transcritical (pitchfork) bifurcation happens with points A, B , 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 , 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, ). 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 . 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 . 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?
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. ) and temporal  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 , 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 . 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.
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.
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.
We thank Alex Lobkovsky for helpful discussions. YIW, EVK and GPK are supported through intramural funds of the US Department of Health and Human Services (to the National Library of Medicine).
- Stern A, Sorek R: The phage-host arms race: shaping the evolution of microbes. Bioessays. 2011, 33 (1): 43-51.PubMedPubMed CentralView ArticleGoogle Scholar
- Lotka AJ: Elements of Physical Biology. 1925, NY: Williams and WilkinsGoogle Scholar
- Volterra V: Variazioni e fluttuazioni del numero d’individui in specie animali conviventi. Mem Acad Lincei Roma. 1926, 2: 31-113.Google Scholar
- May RM: Stability and Complexity in Model Ecosystems. 2001, Princeton: Princeton University PressGoogle Scholar
- Volterra V: Le cons sur la Theorie Mathematique de la Lutte popur la Vie. 1931, Paris: Gauthier VillareGoogle Scholar
- CRISPR-Cas Systems. RNA-mediated Adaptive Immunity in Bacteria and Archaea. Edited by: Barrangou R, Oost J Van der. 2013, Heidelberg: SpringerGoogle Scholar
- Wiedenheft B, Sternberg SH, Doudna JA: RNA-guided genetic silencing systems in bacteria and archaea. Nature. 2012, 482 (7385): 331-338.PubMedView ArticleGoogle Scholar
- Marraffini LA, Sontheimer EJ: CRISPR interference: RNA-directed adaptive immunity in bacteria and archaea. Nat Rev Genet. 2010, 11 (3): 181-190.PubMedPubMed CentralView ArticleGoogle Scholar
- Makarova KS, Wolf YI, Snir S, Koonin EV: Defense islands in bacterial and archaeal genomes and prediction of novel defense systems. J Bacteriol. 2011, 193 (21): 6039-6056.PubMedPubMed CentralView ArticleGoogle Scholar
- Koonin EV, Makarova KS: CRISPR-Cas: Evolution of an RNA-based adaptive immunity system in prokaryotes. RNA Biol. 2013, 10 (5): 679-686.PubMedPubMed CentralView ArticleGoogle Scholar
- Sorek R, Lawrence CM, Wiedenheft B: CRISPR-mediated adaptive immune systems in bacteria and archaea. Annu Rev Biochem. 2013, 82: 237-266.PubMedView ArticleGoogle Scholar
- Horvath P, Barrangou R: CRISPR/Cas, the immune system of bacteria and archaea. Science. 2010, 327 (5962): 167-170.PubMedView ArticleGoogle Scholar
- Koonin EV, Wolf YI: Is evolution Darwinian or/and Lamarckian?. Biol Direct. 2009, 4: 42-PubMedPubMed CentralView ArticleGoogle Scholar
- Koonin EV: The Logic of Chance: The Nature and Origin of Biological Evolution. Upper Saddle River. 2011, Upper Saddle River, NJ: FT pressGoogle Scholar
- Sampson TR, Weiss DS: Exploiting CRISPR/Cas systems for biotechnology. Bioessays. 2014, 36 (1): 34-38.PubMedPubMed CentralView ArticleGoogle Scholar
- Ran FA, Hsu PD, Wright J, Agarwala V, Scott DA, Zhang F: Genome engineering using the CRISPR-Cas9 system. Nat Protoc. 2013, 8 (11): 2281-2308.PubMedPubMed CentralView ArticleGoogle Scholar
- Gasiunas G, Siksnys V: RNA-dependent DNA endonuclease Cas9 of the CRISPR system: Holy Grail of genome editing?. Trends Microbiol. 2013, 21 (11): 562-567.PubMedView ArticleGoogle Scholar
- Mali P, Yang L, Esvelt KM, Aach J, Guell M, Dicarlo JE, Norville JE, Church GM: RNA-guided human genome engineering via Cas9. Science. 2013, 339 (6121): 823-826.PubMedPubMed CentralView ArticleGoogle Scholar
- Jiang W, Bikard D, Cox D, Zhang F, Marraffini LA: RNA-guided editing of bacterial genomes using CRISPR-Cas systems. Nat Biotechnol. 2013, 31 (3): 233-239.PubMedPubMed CentralView ArticleGoogle Scholar
- Cho SW, Kim S, Kim JM, Kim JS: Targeted genome engineering in human cells with the Cas9 RNA-guided endonuclease. Nat Biotechnol. 2013, 31 (3): 230-232.PubMedView ArticleGoogle Scholar
- Fitch WM: The variety of human virus evolution. Mol Phylogenet Evol. 1996, 5 (1): 247-258.PubMedView ArticleGoogle Scholar
- Campbell AM: Conditions for the existence of bacteriophage. Evolution. 1960, 15: 153-165.View ArticleGoogle Scholar
- Levin BR, Stewart FM, Chao L: Resource-limited growth, competition, and predation: A model and experimental studies with bacteria and bacteriophage. Am Nat. 1977, 111: 3-24.View ArticleGoogle Scholar
- Williams HT: Phage-induced diversification improves host evolvability. BMC Evol Biol. 2013, 13: 17-PubMedPubMed CentralView ArticleGoogle Scholar
- Held NL, Herrera A, Cadillo-Quiroz H, Whitaker RJ: CRISPR associated diversity within a population of Sulfolobus islandicus. PLoS One. 2010, 5 (9): e12988-PubMedPubMed CentralView ArticleGoogle Scholar
- Weinberger AD, Sun CL, Pluciński MM, Denef VJ, Thomas BC, Horvath P, Barrangou R, Gilmore MS, Getz WM, Banfield JF: Persisting Low-abundance viral sequences shape microbial CRISPR-based immunity. PLoS Comp Biol. 2012, in pressGoogle Scholar
- Levin BR, Moineau S, Bushman M, Barrangou R: The population and evolutionary dynamics of phage and bacteria with CRISPR-mediated immunity. PLoS Genet. 2013, 9 (3): e1003312-PubMedPubMed CentralView ArticleGoogle Scholar
- Weinberger AD, Wolf YI, Lobkovsky AE, Gilmore MS, Koonin EV: Viral diversity threshold for adaptive immunity in prokaryotes. MBio. 2012, 3 (6): e00412-e00456.View ArticleGoogle Scholar
- Childs LM, Held NL, Young MJ, Whitaker RJ, Weitz JS: Multiscale model of CRISPR-induced coevolutionary dynamics: diversification at the interface of Lamarck and Darwin. Evolution. 2012, 66 (7): 2015-2029.PubMedPubMed CentralView ArticleGoogle Scholar
- Levin BR: Nasty viruses, costly plasmids, population dynamics, and the conditions for establishing and maintaining CRISPR-mediated adaptive immunity in bacteria. PLoS Genet. 2010, 6 (10): e1001171-PubMedPubMed CentralView ArticleGoogle Scholar
- Gandon S, Vale PF: The evolution of resistance against good and bad infections. J Evol Biol. 2014, 27 (2): 303-312.PubMedView ArticleGoogle Scholar
- Iranzo J, Lobkovsky AE, Wolf YI, Koonin EV: Evolutionary dynamics of the prokaryotic adaptive immunity system CRISPR-Cas in an explicit ecological context. J Bacteriol. 2013, 195 (17): 3834-3844.PubMedPubMed CentralView ArticleGoogle Scholar
- Han P, Niestemski LR, Barrick JE, Deem MW: Physical model of the immune response of bacteria against bacteriophage through the adaptive CRISPR-Cas immune system. Phys Biol. 2013, 10 (2): 025004-PubMedView ArticleGoogle Scholar
- Haerter JO, Sneppen K: Spatial structure and Lamarckian adaptation explain extreme genetic diversity at CRISPR locus. MBio. 2012, 3 (4): e00112-e00126.View ArticleGoogle Scholar
- Haerter JO, Trusina A, Sneppen K: Targeted bacterial immunity buffers phage diversity. J Virol. 2011, 85 (20): 10554-10560.PubMedPubMed CentralView ArticleGoogle Scholar
- Bikard D, Marraffini LA: Innate and adaptive immunity in bacteria: mechanisms of programmed genetic variation to fight bacteriophages. Curr Opin Immunol. 2012, 24 (1): 15-20.PubMedView ArticleGoogle Scholar
- Rechavi O: Guest list or black list: heritable small RNAs as immunogenic memories. Trends Cell Biol. 2013, 24 (4): 212-220.PubMedView ArticleGoogle Scholar
- Rimer J, Cohen IR, Friedman N: Do all creatures possess an acquired immune system of some sort?. Bioessays. 2014, 36 (3): 273-281.PubMedView ArticleGoogle Scholar
- Fineran PC, Charpentier E: Memory of viral infections by CRISPR-Cas adaptive immune systems: acquisition of new information. Virology. 2012, 434 (2): 202-209.PubMedView ArticleGoogle Scholar
- Westra ER, Swarts DC, Staals RH, Jore MM, Brouns SJ, van der Oost J: The CRISPRs, they are a-changin’: how prokaryotes generate adaptive immunity. Annu Rev Genet. 2012, 46: 311-339.PubMedView ArticleGoogle Scholar
- Kuznetsov YA: Elements of applied bifurcation theory. 1995, Vienna: SpringerView ArticleGoogle Scholar
- Makarova KS, Wolf YI, Koonin EV: Comparative genomics of defense systems in archaea and bacteria. Nucleic Acids Res. 2013, 41 (8): 360-4377.View ArticleGoogle Scholar
- Makarova KS, Wolf YI, Koonin EV: The basic building blocks and evolution of CRISPR-cas systems. Biochem Soc Trans. 2013, 41 (6): 1392-1400.PubMedView ArticleGoogle Scholar
- Aviram I, Rabinovitch A: Bifurcation analysis of bacteria and bacteriophage coexistence in the presence of bacterial debris. Commun Nonlinear Sci Numer Simul. 2012, 17: 242-254.View ArticleGoogle Scholar
- Aviram I, Rabinovitch A: Bacteria and lytic phage coexistence in a chemostat with periodic nutrient supply. Bull Math Biol. 2014, 76 (1): 225-244.PubMedView ArticleGoogle Scholar
- Juliano C, Wang J, Lin H: Uniting germline and stem cells: the function of Piwi proteins and the piRNA pathway in diverse organisms. Annu Rev Genet. 2011, 45: 447-469.PubMedView ArticleGoogle Scholar
- Stuwe E, Toth KF, Aravin AA: Small but sturdy: small RNAs in cellular memory and epigenetics. Genes Dev. 2014, 28 (5): 423-431.PubMedPubMed CentralView ArticleGoogle Scholar
- Landry CD, Kandel ER, Rajasethupathy P: New mechanisms in memory storage: piRNAs and epigenetics. Trends Neurosci. 2013, 36 (9): 535-542.PubMedView ArticleGoogle Scholar
- Luteijn MJ, Ketting RF: PIWI-interacting RNAs: from generation to transgenerational epigenetics. Nat Rev Genet. 2013, 14 (8): 523-534.PubMedView ArticleGoogle Scholar
- Paez-Espino D, Morovic W, Sun CL, Thomas BC, Ueda K, Stahl B, Barrangou R, Banfield JF: Strong bias in the bacterial CRISPR elements that confer immunity to phage. Nat Commun. 2013, 4: 1430-PubMedView ArticleGoogle Scholar
- Delaney NF, Balenger S, Bonneaud C, Marx CJ, Hill GE, Ferguson-Noel N, Tsai P, Rodrigo A, Edwards SV: Ultrafast evolution and loss of CRISPRs following a host shift in a novel wildlife pathogen, Mycoplasma gallisepticum. PLoS Genet. 2012, 8 (2): e1002511-PubMedPubMed CentralView ArticleGoogle Scholar
- Kuno S, Sako Y, Yoshida T: Diversification of CRISPR within coexisting genotypes in a natural population of the bloom-forming cyanobacterium Microcystis aeruginosa. Microbiology. 2014, 160 (Pt 5): 903-916.PubMedView ArticleGoogle Scholar
- Makarova KS, Haft DH, Barrangou R, Brouns SJ, Charpentier E, Horvath P, Moineau S, Mojica FJ, Wolf YI, Yakunin AF, van der Oost J, Koonin EV: Evolution and classification of the CRISPR-Cas systems. Nat Rev Microbiol. 2011, 9 (6): 467-477.PubMedView ArticleGoogle Scholar
- Bazykin AD: Nonlinear Dynamics of Interacting Populations, Volume 11. 1998, Singapore, New Jersey, London, Hong Kong: World ScientificGoogle Scholar
- Khibnik AI, Kuznetsov YA, Levitin VV, Nikolaev EV: Continuation techniques and interactive software for bifurcation analysis of ODEs and iterated maps. Physica D. 1993, 62: 360-371.View ArticleGoogle Scholar
- Heilmann S, Sneppen K, Krishna S: Coexistence of phage and bacteria on the boundary of self-organized refuges. Proc Natl Acad Sci U S A. 2012, 109 (31): 12828-12833.PubMedPubMed CentralView ArticleGoogle Scholar
- Maslov S, Sneppen K: Well-temperate phage. arXiv:13081646 [q-bioPE]. 2013Google Scholar
- Weitz JS, Dushoff J: Alternative stable states in host-phage dynamics. Theor Ecol. 2008, 1: 13-19.View ArticleGoogle Scholar
- Wang Z, Goldenfeld N: Fixed points and limit cycles in the population dynamics of lysogenic viruses and their hosts. Phys Rev E Stat Nonlin Soft Matter Phys. 2010, 82 (1 Pt 1): 011918-PubMedView ArticleGoogle Scholar
- Edgar R, Qimron U: The Escherichia coli CRISPR system protects from lambda lysogenization, lysogens, and prophage induction. J Bacteriol. 2010, 192 (23): 6291-6294.PubMedPubMed CentralView ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.