We model bacteria and phage interacting and coevolving in a well-mixed system (Figure 1A and ‘Model’). Bacteria divide by consuming nutrients and phages reproduce by creating a burst of $B$ new phages after successfully infecting a bacterium. Bacteria can contain a single CRISPR spacer that confers immunity against phages with a matching protospacer. Phages are labelled with a single protospacer type, a binary sequence of length $L=30$ that can mutate to a new type during a burst with probability $\mu L$, where $\mu$ is the per-base mutation rate per generation. All simulations begin with a single clonal phage population unless otherwise specified.

Figure 1 with 7 supplements see all

Download asset Open asset

Coexistence occurs across a wide range of parameters but is not guaranteed: below a certain success probability $p_V^0=\frac{1}{B}\left(\frac{gf}{(1-f)\alpha }+1\right)$, phages are not able to reproduce often enough to overcome their base death rate due to outflow and adsorption and are driven extinct (grey area of Figure 1D; Bonsma-Fisher et al., 2018). In this expression, $g$ is the bacterial growth rate, $B$ is the phage burst size, $\alpha$ is the phage adsorption rate, $f=F/(g{C}_0)$ is a normalized outflow rate, and C₀ is the inflow nutrient concentration. This is the same extinction threshold reported by Payne et al., 2018 as the cutoff for achieving herd immunity in a well-mixed bacterial population. To a first approximation, phages must successfully infect every $1/B$ bacteria they encounter, but if bacteria are growing quickly, then phages must do better to overcome bacterial growth, leading to the extra terms in this expression (see ‘Phage extinction threshold’). We write this extinction threshold as $A=\frac{(B{p}_V-1)(1-f)\alpha }{fg}$. Above the phage extinction threshold (${\displaystyle A>1}$), the phage population size increases with increasing $p_V$ but eventually decreases again as bacterial numbers are driven too low to support a large phage population (Bonsma-Fisher et al., 2018). A similar non-monotonicity as a function of the probability of naive bacterial resistance ($1-p_V$) was described in theoretical work by Weinberger et al., 2012b. In our model the position of the peak in phage population size as a function of infection success probability is determined by $e$, the effectiveness of CRISPR spacers against phage; increasing $e$ pushes the peak to higher $p_V$ (Figure 1D). While $e$ is a constant parameter that determines the outcome of pairwise interactions between bacteria and phages, the bacterial population as a whole possesses an average immunity to phages that is a weighted average of all the possible pairwise interactions (Figure 1A inset). It is the overall average immunity that determines population outcomes, which we describe in detail in ‘Pathogen and host diversity must be considered together’.

To focus on regimes where bacteria and phages coexist, we select parameters within the deterministic coexistence regime to explore bacteria-phage coevolution. Even in this regime, stochastic extinction will eventually come for one or both populations in simulations (Figure 1E), though the timescale of extinction may be extremely long for large population sizes (Badali and Zilman, 2020). Phages are more susceptible to stochastic extinction than bacteria because of their large burst size $B$ which increases their overall population fluctuations (Appendix 3). The length of coexistence before stochastic extinction depends on population size as well as simulation parameters and initial conditions (see ‘Stochastic population extinction’): phage populations can be rescued from extinction by high mutation rate (Figure 1—figure supplement 3) or high initial protospacer diversity (Figure 1—figure supplement 6), but are more likely to go extinct if spacer effectiveness is high (Figure 1—figure supplement 1). Conversely, bacteria are more likely to go extinct if spacer effectiveness or spacer acquisition rate are low (Figure 1—figure supplement 4). Population survival and persistence in natural populations is impacted by additional factors we do not address in our model, including immigration (Volkov et al., 2003; Chabas et al., 2016, niche partitioning Simek et al., 2010; Weitz et al., 2013; Mills et al., 2013; Badali and Zilman, 2020; Voigt et al., 2021), environmental fluctuations (Abreu et al., 2020; Voigt et al., 2021), and spatial structure (Haerter et al., 2011; Haerter and Sneppen, 2012; Heilmann et al., 2010; Heilmann et al., 2012; Simmons et al., 2018; Skanata and Kussell, 2021).

Across a wide range of coexistence parameters, our simulations show continual phage evolution and bacterial CRISPR adaptation in response (Figure 1B). New phage protospacer clones arise from a single founding clone by mutation, and a small fraction of new mutants grow to a large size and become established. Once phage clones become large, bacteria acquire matching spacers and an immune bacterial subpopulation becomes established. The specific protospacer and spacer types present in the population continually change as old types go extinct and new types are created by phage mutation, but the average total diversity and average overlap between bacteria and phage remains constant at steady state (Figure 8). Both bacteria and phage clones stochastically go extinct, completing the life cycle of a clonal population (Figure 1C).

New phage protospacer clones continually arise and go extinct in our simulations, generating turnover in clone identity in the population. Despite constant turnover, however, the total number of clones remains fixed at steady state. We use the mean number of bacterial clones at steady state, designated $m$, as a measure of system diversity. This choice of diversity measurement is equivalent to the Hartley entropy of the clone size distribution, a special case of the Rényi entropy (Mora and Walczak, 2016; Altan-Bonnet et al., 2020). This definition weights all clones equally regardless of their abundance; such a measurement is not appropriate when clone size distributions are very broad and small clones may be unsampled, but is reasonable when clone size distributions are relatively narrow and all clones are sampled (Mora and Walczak, 2016). In our simulation results, both bacteria and phage populations exhibit relatively narrow clone size distributions across a range of parameters, with the exception of low values of spacer acquisition $\eta$ (Figure 2A, Figure 2—figure supplement 3, Figure 2—figure supplement 4). Even at low $\eta$, however, clone size distributions are approximately exponential, indicating that they are not scale-invariant and that the mean clone size still captures important information about the full clone size distribution.

Figure 2 with 4 supplements see all

Download asset Open asset

What determines clonal diversity? Many factors that correlate with transient diversity have been experimentally identified, such as phage extinction and slower phage evolution at high bacterial spacer diversity (van Houte et al., 2016; Common et al., 2020) and maintenance of a diverse bacterial population when exposed to diverse phages (Paez-Espino et al., 2015; Common et al., 2019; Guillemet et al., 2021; Lopatina et al., 2019), but a conceptual framework to understand emergent diversity has remained elusive. For instance, while initial high spacer diversity puts low-diversity phage populations under intense pressure to the point of driving them extinct (van Houte et al., 2016; Common et al., 2020), is the same true for emergent bacterial diversity after an extended period of coexistence? Is observed high bacterial spacer diversity indicative of successful bacterial escape from phage predation or an indicator of increased phage pressure? In our model, phage and bacterial diversity is tightly coupled: the number of large phage clones is approximately the same as the number of bacterial clones (Figure 22). This is also the case in experimental coevolution data from Paez-Espino et al., 2015: the number of phage protospacer types is on the same order of magnitude as the number of bacterial spacer types across most similarity thresholds (Figure 72). There is evidence that this coupling of diversity may also occur in the wild: a recent longitudinal study of Gordonia bacteria interacting with phage in a wastewater treatment plant identified 14 high-coverage phage genotypes and 11 high-coverage bacterial variants based on CRISPR spacer sequence (Guerrero et al., 2021a).

Using the tight correspondence between bacterial diversity and phage diversity in our model, we calculate the overall steady-state diversity by balancing the effective phage clone mutation rate $\overline{\mu }=\frac{1}{g{C}_0}\alpha B{p}_V(1-{\text{e}}^{-\mu L})(1-\frac{e\nu }{m})n_V{n}_B$, phage clone establishment probability $P_{est}$, and the time to extinction for large phage clones $T_{ext}$ (details in ‘Measuring diversity):

Equation 1 arises from the simple statement that the number of large clones must be equal to their establishment rate ($P_{est}\overline{\mu }$) multiplied by their average time to extinction ($T_{ext}$). This relationship successfully predicts the number of bacterial clones at steady state across a wide range of parameters and a wide range of diversity values (Figure 2B inset, Figure 2—figure supplement 1, and Figure 2—figure supplement 2). At the lowest value of $\eta$ the prediction tends to overestimate the number of clones – in this regime, low acquisition means that phage clones go extinct because of clonal interference before bacteria are able to acquire spacers (Figure 16, ‘Measuring diversity’).

Through approximations (‘Analytic approximations for diversity’), we find that diversity depends on a single combined parameter to the power $1/3$ (Figure 2B, Equation 2), and this parameter is proportional to spacer effectiveness $e$, the probability of bacterial survival followed by spacer acquisition $(1-p_V)\eta$, and the phage mutation rate $\mu L$.

Each of these parameters intuitively increases diversity (e.g., a higher phage mutation rate means that phage diversity increases and bacterial diversity follows suit). What is surprising is that their combined effect on diversity is to a power much less than 1: this 1/3 exponent means that if the mutation rate increased tenfold, the diversity would only increase by about a factor of two (Appendix 2—figure 1). In contrast, a simple neutral model of cell division with mutations gives a linearly proportional increase in diversity for the same increase in mutation rate (Appendix 2—figure 2).

To understand where the dependence of diversity on these key parameters come from, we look more closely at each component expression. The effective phage mutation rate ${\displaystyle \overline{\mu }}$ depends linearly on the parameter ${\displaystyle \mu :\overline{\mu }\approx \frac{g{C}_0(1-f)f\mu L}{\alpha P_V}}$, while both the probability of establishment and the time to extinction depend inversely on diversity ${\displaystyle m:{\mathit{p}}_{est}\approx \frac{1}{m}\frac{ce\eta (1-P_V)g{C}_0(1-f)}{B_{p_{V}r}}}$ (Equation 5) and ${\displaystyle T_{ext}\approx \frac{1}{m}\frac{2g{C}_0(1-f)}{f\alpha B{p}_V}}$ (Equation 173, Appendix 3). By comparison with Equation 1, we find that $m^3$ depends approximately linearly on mutation rate, resulting in the weak $m\propto {\mu }^{\frac{1}{3}}$ dependence on mutation rate.

The dependence of diversity on both $e$ and $\eta$ comes from the probability of phage establishment since $\overline{\mu }$ depends only very weakly on these parameters through its dependence on total population sizes and $T_{ext}$ depends explicitly on $m$ alone, not $e$ or $\eta$. The phage probability of establishment is proportional to $\frac{e\eta (1-p_V)}{m}$ (Equation 5), and as before, this gives $m^3\propto e\eta (1-p_V)$. Bacteria are more successful at high $\eta (1-p_V)e$, which increases the phage establishment probability. Previous theoretical work has predicted that diversity increases as spacer acquisition rate increases (Childs et al., 2012); here, we provide a quantitative prediction for this dependence. In the following sections, we explore phage establishment in more detail.

We find that diversity emerges in our model from the balance of phage clone establishment and extinction. However, only some phage mutants escape initial stochastic extinction and survive long enough to become established. What determines the fate of a new phage mutant? In our model, a single phage mutation in a protospacer can completely overcome CRISPR targeting, which means that new phage mutants can infect all bacteria equally well and their initial growth rate s₀ is independent of CRISPR: $s_0\approx \alpha n_B(B{p}_V-1)-F$, where $F$ is the chemostat flow rate (a shared death rate for phages and bacteria). Surprisingly, however, even once bacteria start to acquire matching spacers, the probability of establishment for new phage mutants is still well-described by theory in which CRISPR targeting only influences total average population sizes (Figure 3C); that is, the specific interaction between a phage and its matching clone can be ignored. Intuitively, this is because phage clones must grow to a certain size before bacteria encounter them enough to begin to acquire spacers, and this size turns out to be large enough to avoid stochastic extinction (Figure 15). The probability of phage establishment is $\frac{2s_0}{B(s_0+{\delta }_0)}$, where ${\delta }_0=F+\alpha n_B(1-p_V)$ is the initial phage mutant death rate. Importantly, these rates are independent of the population size of matching CRISPR bacteria clones; the only dependence on bacteria is on the total bacterial population size $n_B$ (Appendix 3—figure 1), which is fairly stable at steady state.

Figure 3 with 3 supplements see all

Download asset Open asset

Even though CRISPR targeting does not explicitly affect the establishment probability for new phage mutants, the probability of phage establishment increases as the average bacterial immunity increases (Figure 3C). Average immunity is a measure of the overall effectiveness of CRISPR immunity for the entire bacterial population; it is the average of all pairwise immunities between phage clones and bacterial clones weighted by their population sizes (Figure 1A inset). Because higher average immunity is beneficial for bacteria and leads to larger bacterial population sizes (Figure 3C), higher average immunity also means there is stronger selective advantage for new phage mutants that can escape CRISPR targeting (Figure 3C and Figure 5A).

For insight into the nature of this dependence, we approximated the probability of establishment at two extremes of average immunity (see ‘Dominant balance approximations for $\nu$’). The probability of establishment can also be written

where $\nu$ is the fraction of bacteria that have spacers and $m$ is the average number of bacterial clones at steady state. Low average immunity occurs when $e/m$ is small: in this limit, the fraction of bacteria with spacers is

where $A=\frac{(B{p}_V-1)(1-f)\alpha }{fg}$ is the extinction threshold for phages (${\displaystyle A>1}$ for phage survival), and ${\tilde{n}}_V$ is the phage population size calculated with $e=0$ (Equation 52); this is the extreme limit of low average immunity where total population sizes approach their values in the absence of CRISPR immunity (Figure 26). The quantity $\frac{r}{\eta (1-p_V)\alpha {\tilde{n}}_V}$ can be understood as the ratio of the rate of spacer loss to the rate of spacer acquisition. If spacer loss is high and acquisition is low, the fraction of bacteria with spacers ($\nu$) decreases. This in turn correlates with a lower probability of phage establishment (Equation 3) since phages are less likely to establish when there is less pressure from bacteria defending against them with CRISPR. Equation 4 is plotted in Figure 3C in blue. At very high average immunity, on the other hand, $\nu$ becomes close to 1 regardless of $\eta$ (Figure 3—figure supplement 2), which is why the curves in Figure 3C overlap at high average immunity ($\nu =1$ is plotted in green).

For more insight into the parameter dependence of $P_{est}$, we apply the same approximation as for $m$, expanding $P_{est}$ in $1/B$ and $\eta$ (‘Approximation for $m$’):

The probability of establishment increases with the probability of bacterial escape and spacer acquisition $\eta (1-p_V)$ and with spacer effectiveness $e$, but decreases with increasing mutation rate $\mu L$. This is consistent with the intuition that more successful bacteria increase the strength of selection for phage mutants (higher $\eta$ and $e$), but that higher mutation rate reduces the probability of establishment for any particular mutant.

We next asked how cross-reactivity between spacer and protospacer types impacts population dynamics and outcomes. Adding cross-reactivity was motivated by several experimental observations in CRISPR immunity: (a) In type I and type II CRISPR systems, single mutations in the PAM or protospacer seed regions (approximately 8 nucleotides at the start of a protospacer) can facilitate phage escape, whereas mutations elsewhere in the protospacer are tolerated by the CRISPR system (Deveau et al., 2008; Pyenson and Marraffini, 2020). Even when a phage manages to escape direct targeting, in type I and II systems an imperfect spacer match can facilitate priming: when Cas machinery binds to a protospacer match, even if unable to cleave the target, the likelihood of acquiring a nearby spacer is increased (Westra et al., 2015; Rao et al., 2017; Pyenson and Marraffini, 2020; Weissman et al., 2020). (b) Type III CRISPR targeting in Staphylococcus epidermidis has been shown to be completely tolerant to all single and even double mutations, meaning that these systems are naturally cross-reactive (Pyenson et al., 2017).

We simulated two types of cross-reactivity: one in which phages experience an exponential decrease in CRISPR effectiveness with mutational distance (Equation 7), and a step-function cross-reactivity where phages require additional $1\le \theta \le 3$ mutations to perfectly escape CRISPR targeting (Equation 8). Phage success probability without cross-reactivity is given by Equation 6, a special case of 7 and 8 with $\theta =0$.

The mutational distance $d_{ij}=\sum |(V_i-V_j)|$ is the number of mutations between a protospacer and spacer, and $\theta$ scales the radius of cross-reactivity, with larger $\theta$ meaning more mutations are required to achieve the same immune escape (Yan et al., 2019).

Exponential cross-reactivity has been modelled extensively in vertebrate adaptive immunity (Chao et al., 2005; Mayer et al., 2015; Marchi et al., 2019; Yan et al., 2019; Schnaack and Nourmohammad, 2021; Chardès et al., 2022), and our definition of exponential cross-reactivity as a function of mutational distance is the same as in Yan et al., 2019. Step-function cross-reactivity, on the other hand, is reminiscent of the type III CRISPR system in which multiple point mutations are required to escape CRISPR targeting (Pyenson et al., 2017) and was also modelled theoretically by Han et al., 2013.

We find that cross-reactivity results in strikingly different dynamics of clone establishment and persistence (Figure 4 and ‘Cross-reactivity’), including a travelling wave regime in which genetically neighbouring clones ‘pull’ each other along (Figure 4B–D), giving way to a regime in which new mutants have very low establishment rates in the case of step-function cross-reactivity (Figure 4C). When cross-reactivity is high and phages require multiple mutations to escape targeting, we expect new phage mutants to have a lower fitness on average because they may already be within the immunity range of existing bacterial clones. However, unlike without cross-reactivity, not all new mutant fitnesses are the same because fitness now depends on the distribution of matching spacers in the bacterial population. Cross-reactivity adds a ‘direction’ in the genetic landscape and a fitness gradient for new phage mutants, leading to a series of rapid establishments and a travelling-wave regime. We believe these rapid establishments also suppress instantaneous diversity; diversity emerges longitudinally instead of concurrently in this regime (Figure 4 and Figure 61). This is true for both types of cross-reactivity we studied, but we also observed striking qualitative differences between exponential and step-function cross-reactivity. In the former, all phage mutants are guaranteed a slight fitness advantage, and the travelling wave appears right at the start of simulations, with occasional lineage-splitting events (Figure 4B). In contrast, step-function cross-reactivity means that mutants within the cross-reactivity radius have no fitness advantage at all. In this case, there is some variable initial length of time required to establish the travelling wave pattern: at the start of a simulation, all mutants are within the cross-reactivity radius and evolve purely neutrally. At least $d=\theta$ new mutants must establish through neutral dynamics before the next mutant can escape CRISPR targeting and grow with high fitness; this leads to the travelling-wave regime appearing later for $\theta =2$ than for $\theta =1$, and sometimes never appearing at all. Once the traveling wave appears, the dominant phage and bacteria clone types are exactly offset by $\theta$: the most abundant phage clone will in general be $\theta$ mutations away from the most abundant concurrent bacteria clone. This is directly visible in the clone trajectories in Figure 4C–D inset and Figure 54: matching-colour bacteria and phage clones are offset in time. Another regime is also possible in the case of step-function cross-reactivity: if multiple large clones establish that are each outside of all the others’ cross-reactivity radii (Figure 63), then new mutants all have zero fitness and the travelling wave comes to an abrupt halt (around 7000 generations in Figure 4C). Establishment of new clones is once again extremely rare, and the existing large clones persist for a long time (Figures 53 and 62).

Figure 4 with 4 supplements see all

Download asset Open asset

Changes to the fitness landscape of individual clones caused by cross-reactivity also influence population-level outcomes such as diversity. In general, we find that for high cross-reactivity, average diversity is lower than predicted for simulations without cross-reactivity because the probability of phage establishment per mutant decreases (Figure 5A). A decrease in phage mutant fitness as the strength of cross-reactivity increases is in line with the dependence of fitness on cross-reactivity reported by Yan et al., 2019 and Rouzine and Rozhnova, 2018. Even with more complicated underlying dynamics, though, measuring average immunity alone is enough to reproduce total population sizes using our simple deterministic equations developed in the limit of no cross-reactivity (Figure 5C, 'Total population size'). This is because average immunity completely captures the population-level impact of CRISPR immunity, and we can replace the CRISPR effectiveness parameter $e$ with measured average immunity in Equations 13–17 to get very good agreement between measured and predicted population sizes even away from steady state (Figure 5—figure supplement 1). In contrast, inferring average immunity from the number of clones using the approximation that effective $e\approx e/m$ does not give good agreement with simulation results for simulations with cross-reactivity (Figure 5B).

Figure 5 with 1 supplement see all

Download asset Open asset

Figure 5B makes an additional subtle point: higher immune diversity promotes larger population sizes for phages but not for bacteria. This is counterintuitive and appears to conflict with several experimental results that show that more bacterial diversity increases the likelihood of phage extinction (van Houte et al., 2016; Common et al., 2020) and decreases the ability of phages to adapt (Morley et al., 2017). This discrepancy appears at first glance to be a consequence of limiting bacteria and phages to a single spacer and protospacer in our model, but we propose that it is actually because bacteria and phage diversity is decoupled in prior experiments. In the following section, we introduce a toy model to understand the conceptual impact of multiple protospacers and spacers when overall phage and bacterial diversity is coupled.

Our model restricts bacteria to a single spacer and phages to a single protospacer. This single-spacer assumption is a reasonable approximation for the experimental systems we consider: many liquid culture experiments find that most bacteria acquire a single spacer when exposed to phages, (Heler et al., 2015; Heler et al., 2019; Pyenson and Marraffini, 2020) even after up to 2 weeks (Paez-Espino et al., 2013; Common et al., 2019). Additionally, metagenomic results show that most spacer matches to viral sequences are located near the leader end of the CRISPR array Weinberger et al., 2012a. These results, combined with theory positing that recently acquired spacer provide more immune benefit than older spacers (Childs et al., 2012; Han et al., 2013), suggest that dynamics may be dominated by a very small number of spacers per bacterium.

Nevertheless, many bacteria possess tens to hundreds of spacers (Pavlova et al., 2021) and the question of multiple protospacers remains. Can we make inferences about how a more realistic multiple spacer or multiple protospacer scenario will change the relationship between diversity and average immunity in our model? First, we note that our model with exponential cross-reactivity is a good approximation for a realistic scenario where bacteria are limited to a single spacer but phages can have multiple protospacers. In this situation, multiple bacteria with unique spacers may target different protospacers in the same phage. A mutation in one protospacer will increase phage fitness but will not provide perfect immune escape, just as in the case of exponential cross-reactivity (Figure 4B). Even in this scenario, our qualitative results remain the same: bacteria and phage diversity is coupled and bacteria have higher average immunity at low diversity (Figure 5B and Figure 50).

We explored the consequences of both multiple spacers and multiple protospacers with a toy model. We generated synthetic sets of protospacers and spacers and varied the total diversity by changing the size of their pool and distributing their abundances exponentially to qualitatively match data (Bonsma-Fisher et al., 2018). We randomly assigned protospacers and spacers to individual phages and bacteria, then calculated average immunity. Importantly, we constrain phage and bacterial diversity to be coupled – overall spacer and protospacer diversity is the kept the same in all scenarios.

We find that across all combinations of array sizes and for two of three scenarios of protospacer diversity, average immunity is negatively correlated with diversity (Figure 69). Only when we limited phage diversity to a single mutating protospacer with all other protospacers conserved were bacteria with multiple spacers reliably able to target conserved protospacers regardless of the level of diversity in the variable protospacer position (Figure 69D). We compared our toy model to the experimental setup from Common et al., 2020. In their experiments, bacterial CRISPR spacer diversity was increased while phage diversity was kept constant with a single phage strain able to infect only one of the bacterial clones. This translated to an increasing initial average immunity as bacterial diversity increased and phages were able to infect a smaller proportion of the bacterial population (Figure 69A), which we think is not a realistic situation in a coevolving population of phages and bacteria.

Based on these results, we predict that in many realistic situations average immunity does not increase with diversity if phage and bacterial diversity are coupled. The intuition is this: if diversity of both bacteria and phages increase beyond the functional length of the CRISPR array, bacteria cannot be immune to all phage strains at once and average immunity must go down as diversity increases. The actual benefit of bacterial diversity depends very strongly on the characteristics of phage diversity (Figure 69) and manipulations of diversity should be understood in terms of their impact on average immunity. This has implications for understanding the so-called 'dilution effect' which is typically described as a decrease in the fraction of susceptible hosts as host diversity increases (Chabas et al., 2018; Common et al., 2020). Our model and analysis show that a dilution effect can only occur if host diversity increases out of proportion to pathogen diversity as in Common et al., 2020 (Figure 69A), whereas if both host and pathogen diversity increase in tandem, the dilution effect actually changes direction, decreasing the fraction of pathogens that a host may be immune to as pathogen diversity increases (Figure 69B–D). In the same vein, previous theoretical work has found that CRISPR with a fitness cost would be selected against when viral diversity is high for this reason (Weinberger et al., 2012b; Iranzo et al., 2013), notably with models that allowed for multiple spacers and protospacers. Our results, building on existing theory, show that phage and bacterial diversity must be explored in tandem and that CRISPR immunity may provide less benefit when phage diversity is high.

How quickly does the phage population evolve? With explicitly modelled spacer and protospacer sequences of length $L=30$, we quantify the speed of evolution as the ‘mutational distance’ per generation (Rouzine and Rozhnova, 2018): between two time points, how far in genome space has the phage population travelled, or how many mutations have occurred on average (Speed of evolution)?

The speed of evolution is both stable and repeatable between simulations at steady state and highly correlated between bacteria and phage (Figure 38). We find that the speed of evolution is inversely proportional to the time to extinction for large phage clones (Figure 6B, 'Measuring speed of evolution'). Intuitively, if phage clones turn over more quickly (small time to extinction), the population is able to move more quickly to a different genetic state, while if the time to extinction is large, new mutants are seeded close to parent populations that persist for a long time, limiting the mutational distance. Since the phage time to extinction has a simple relationship to bacterial diversity, we can also relate speed to diversity, and we find that the speed of evolution is proportional to diversity and proportional to the same parameter combination to the power 1/3 (Figure 6C). As in the case of diversity, increasing phage mutation rate $\mu$ increases the speed of evolution sublinearly, and like diversity, the prediction diverges from simulation results at low $\eta$ as described in 'Measuring diversity.

Figure 6 with 5 supplements see all

Download asset Open asset

Our simulations show continual spacer turnover at steady state (Figure 6D, Figure 6—figure supplement 3, Figure 6—figure supplement 4), a feature that we might also expect to see in actively evolving laboratory or natural populations of bacteria and phages. We analysed data from a long-term in vitro coevolution experiment with Streptococcus thermophilus bacteria and phage (Paez-Espino et al., 2015), data from a time-series sampling of a natural aquifer community of bacteria (Burstein et al., 2016), and data from a time-series sampling of a wastewater treatment plant (Guerrero et al., 2021a) (see 'Materials and methods') and calculated spacer turnover over time (Figure 6E–G). We found that spacer sequences experienced turnover in all cases, indicating ongoing change in the spacer content in these populations, though no experimental system sho wed complete turnover of all spacer types. Even if spacer loss is happening for non-selective reasons (e.g., mixing or flow in the aquifer system), turnover indicates that new spacers are being acquired as well and that CRISPR systems are active. Moreover, the timescale of early spacer turnover in the S. thermophilus experiment is similar to the timescale in our simulations — after 1000 generations, most simulations have between 40 and 60% of clones remaining (Figure 6D) and about 60% of clones remain after 1000 generations in the experimental data (Figure 6E). We modelled all simulation parameters on known parameters for S. thermophilus where possible.

Bacteria acquire spacers in response to phage clones becoming large. The spacer composition of the bacterial population tracks the phage protospacer composition with a lag: Figure 6A shows the first two components of a PCA decomposition of bacteria and phage abundances over time in a simulation, a visual illustration of bacterial tracking. Without cross-reactivity, trajectories in this lower-dimensional space do not travel in a straight line; they are reminiscent of the diffusive coevolutionary trajectories in antigenic space described in a theoretical model of vertebrate virus-host coevolution (Marchi et al., 2019). In the travelling wave regime with cross-reactivity, trajectories are much more ballistic (Figure 6—figure supplement 1). Can we quantify how well and how quickly can bacteria track the evolving phage population with CRISPR?

Average immunity is a simple metric that quantifies the overlap between bacteria and phages. In experiments with both microbes and vertebrates, time-shift infectivity analyses between host and pathogen populations typically show that hosts are more immune to pathogens from the past and less immune to pathogens from the future (Richman et al., 2003; Frost et al., 2005; Moore et al., 2009; Hall et al., 2011; Koskella, 2014; Betts et al., 2018; Common et al., 2019; Dewald-Wang et al., 2022). We conduct a time-shift analysis on our simulation data (Figure 7A and B, Figure 7—figure supplement 11, Figure 7—figure supplement 12) and find that the same pattern holds true, but only for a limited time window: bacteria are indeed more immune to phages from the past, but this past immunity has a peak and then decays as we look further into the past, eventually reaching 0 immunity (Figure 7B). The presence of a peak in past immunity reflects the timescale of spacer turnover: once the bacterial population has lost all spacers from a previous time point, it is no longer immune to contemporaneous phages and time-shifted average immunity falls to 0. In simulations, the position of this peak is dependent on parameters such as $\eta$ and $\mu$ (Figure 7E and F). As $\mu$ increases, the peak occurs further in the past since phage are now moving more quickly away from bacteria tracking, while as $\eta$ increases, the peak moves closer to the present since bacteria are responding more quickly to changes in the phage population. This suggests that there is a tradeoff between immune memory durability and responding quickly to immune threats: bacteria must choose between tracking the phage population closely or keeping past immunity for a long time.

Figure 7 with 12 supplements see all

Download asset Open asset

Time-shift analyses are usually performed explicitly by directly combining stored samples from different time points (Betts et al., 2018; Common et al., 2019; Dewald-Wang et al., 2022). However, by sequencing CRISPR spacers and phage genomes, a pseudo-analysis can be done without any time-shifted competition experiments by calculating the overlap between bacterial spacers and phage protospacers at different times delays. We performed such an analysis for two published datasets: a long-term laboratory coevolution experiment with S. thermophilus and phage (Paez-Espino et al., 2015) and time series of metagenomic samples over 3 years from a wastewater treatment plant (Guerrero et al., 2021a). In both experiments, whole-genome shotgun sequencing was performed on bacteria and phage DNA, and CRISPR spacers and protospacers were recovered and reported. We re-analysed these datasets to detect CRISPR spacers by finding sequences adjacent to known CRISPR repeats in raw reads. We also detected protospacers by finding matches to our detected CRISPR spacers in reads that did not match the CRISPR repeats or the bacterial reference genome(s) (see ‘Laboratory coevolution experimental data’; ‘Experimental data from wastewater treatment plant’).

We grouped spacers with an 85% similarity threshold and interpolated counts between sequenced time points. We calculated the average overlap assuming $e=1$ (perfect immunity from matching spacers) as a function of time delay, averaging over all combinations of interpolated data with the same time delay (Figure 7C and D). We found two qualitatively different trends as a function of time delay. In the laboratory coevolution experiment, we found that bacteria are more immune to past phages and less immune to future phages (Figure 7C, Figure 7—figure supplement 1), consistent with what we see at short time delay in our model (Figure 7A, black dashed line). In contrast, in the wastewater treatment plant data, we found a peak in average immunity that is roughly centred at zero time delay: bacteria are most immune to phages from their same time point, and immunity rapidly decays in both the past and the future (Figure 7D), qualitatively similar to our simulation results on very long timescales (Figure 7B). These different trends may reflect different regimes of CRISPR immune memory: in the laboratory experiment data, we see no evidence of the decay of immune memory, while in the wastewater treatment data we see a suggestion of decay on the timescale of weeks. Bacteria appear to track the phage population closely in the wastewater data, while in the laboratory data bacteria lag behind the phage population.

Because the variability in average immunity between time points was very high in the wastewater treatment plant data, we performed a Wilcoxon signed-rank test between the average immunity at zero time delay and the average immunity at a time delay ±200 and ±500 (see ‘Calculating average immunity for details’). We found that past immunity after 200 days is significantly lower than present (Wilcoxon statistic $Z=1241$ p-value $p=0.008$), but that immunity after 500 days is not significantly lower than present ($Z=392,p=0.27$). For both time delays, we found that future immunity is significantly lower than present ($Z=580,p=0.0012$ for 500 days and $Z=1293,p=0.003$ for 200 days). Interestingly, these significance values are not symmetric if we pose the question from the perspective of phage: the overlap between phage and future bacteria at 500 days is lower than the overlap for present phages ($Z=507,p=0.024$), and the overlap between phage and past bacteria is lower than the overlap at present ($Z=388,p=0.027$, Figure 88). This asymmetry, that bacteria are generally more immune to all past phages while phage are not more infective against all past bacteria, is qualitatively the same as that reported by Dewald-Wang et al. in an explicit time-shift study of bacteria and phage immunity and infectivity in chestnut trees (Dewald-Wang et al., 2022).

Experimental manipulations of CRISPR diversity have focused on changing the number of bacterial spacers present in the population (van Houte et al., 2016; Morley et al., 2017; Common et al., 2020) while keeping phage diversity fixed (and low), with the exception of Guillemet et al., 2021 who also explored high vs. low phage diversity. In contrast, emergent phage and bacterial diversity are tightly coupled in our model (Figure 2B). At high diversity, bacteria must ‘choose’ which of many phage clones to gain immunity against, meaning that they are then immune to a smaller fraction of the total phage population compared to when diversity is low. This is not a result of limiting bacteria and phage to a single spacer or protospacer; this trend also holds when there are multiple spacers or protospacers provided phage and bacteria diversity are correlated, as we showed in ‘Pathogen and host diversity must be considered together’.

Our toy model described in ‘Pathogen and host diversity must be considered together’ shows that the framework of average immunity provides a conceptually intuitive way to understand the impacts of different types of diversity and modes of evolution. Future work to explicitly model the population dynamics of multiple spacers and protospacers with this lens will allow us to understand how and when these different evolutionary modes arise. For instance, it remains unknown how the linkage of different spacers and protospacers within genomes would affect individual clone dynamics. Our work effectively assumes that all spacers and protospacers evolve independently and are uncoupled, but in real genomes evolving with CRISPR, spacers and protospacers may hitchhike to prominence through selection acting on a different sequence in the same genome. Models of CRISPR locus spacer addition have shown that trailer-end clonality emerges from selective sweeps of newly added spacers (Weinberger et al., 2012a; Han et al., 2013), and hitchhiking will affect interpretations that can be drawn from the abundance of individual spacer clones. Exploring correlations between spacer abundances at different locus positions can be done in experimental data (such as Paez-Espino et al., 2013), and extending our simulation to explicitly include multiple protospacers and spacers per genome will provide further data.

Our results suggest that immune diversity itself may not be the most informative observable feature but that average immunity is what actually determines population outcomes. Indeed, when we introduced cross-reactivity between spacer sequences, we found that diversity decreased and that other processes such as establishment were impacted, but that average immunity still correctly predicted population sizes (Figure 5C). Average immunity completely summarizes the effect of CRISPR immunity and accurately predicts population outcomes. Experiments that manipulate bacterial diversity only change average immunity if they include a changing proportion of sensitive bacteria as in Common et al., 2020 (Figure 69A).

Our results showed complex outcomes and a wide range of overall spacer diversity despite limiting bacteria and phage to a single protospacer and spacer. We also observed rapid loss of previous spacer types as bacteria returned to the naive state before acquiring new spacers. In reality, bacteria can store multiple spacers in their CRISPR arrays (Pavlova et al., 2021; Pourcel et al., 2020) and phages can have several hundred to several thousand protospacer sequences which are determined by a particular protospacer-adjacent motif (PAM) in type I and II CRISPR systems (Leenay et al., 2016). In S. thermophilus phage 2972, for example, there are about 230 protospacers that can be acquired into the CRISPR1 locus of S. thermophilus (Paez-Espino et al., 2013) and 465 that can be acquired into the CRISPR3 locus. Similarly, many bacteria can store tens to hundreds of spacers in their CRISPR arrays (Pavlova et al., 2021; Chen et al., 2022; Bradde et al., 2020; Pourcel et al., 2020). Our single-spacer assumption is reasonable for the experimental systems we consider. In liquid culture experiments with staphylococci, most bacteria acquire a single spacer (Heler et al., 2015; Heler et al., 2019; Pyenson and Marraffini, 2020), and in short-term experiments with S. thermophilus, most bacteria acquired just one CRISPR spacer over 2 weeks (Paez-Espino et al., 2013 or one to three spacers over 9 days Common et al., 2019). In the experimental data we analysed here (a 200+ day experiment), the average number of newly acquired spacers hovered around 5 throughout the experiment (Paez-Espino et al., 2015). In a long-term study of metagenomic data, most spacer matches to extant viral sequences were located near the leader end of CRISPR loci; conserved trailer-end spacers had very few matches to sampled viruses (Weinberger et al., 2012a). All these results suggest that dynamics may be dominated by just a handful of spacers per bacterium. In addition, theoretical work has found that the most recently acquired spacer provides more immune benefit than older spacers (Childs et al., 2012; Han et al., 2013). Our single-protospacer assumption for phages had a more extreme impact, however: in our model, we assumed in the non-cross-reactive case that a single point mutation allowed phages to escape from all CRISPR targeting, when in reality a phage with hundreds of protospacers may need mutations in many different protospacers to fully escape from the bacterial population. Previous work on the impact of CRISPR diversity takes place in this context: that having more diverse bacterial spacers means phages need mutations in multiple protospacers to escape, and that a single mutation only escapes a single bacterial spacer at most, leaving other bacteria still resistant (Chabas et al., 2018; van Houte et al., 2016; Common et al., 2020). We approximated this scenario of gradually increasing phage fitness per protospacer mutation by adding cross-reactivity to our model and showed in ‘Cross-reactivity leads to dynamically unique evolutionary states’ that the relationship between diversity and average immunity does not qualitatively change. We also observed the same qualitative trend in a toy model exploring the effects of both multiple spacers and multiple protospacers (‘Pathogen and host diversity must be considered together’).

Our model makes quantitative predictions about the relationship of population size and population outcomes to average immunity that can be tested in data. We directly calculated average immunity in experimental data without the need to assemble phage genomes (Figure 7), and we found that average immunity is anticorrelated with phage population size in data in one experimental dataset, as predicted by our model (Figure 77). However, we did not find a correlation between average immunity and read count (a proxy for population size) in wastewater treatment plant data (Figure 90). When calculating average immunity using experimental data, one must make assumptions about the immune benefit of spacers. In our work, we assumed that spacer and protospacer sequences that belong to the same similarity group provide perfect immunity. Other more complex assumptions are possible, and there has been a great deal of work quantifying the efficiency of spacers based on the position of SNPs in the protospacer in vitro (see Sternberg et al., 2015 for an example) and in vivo (Soto-Perez et al., 2019). Additionally, phages regularly acquire mutations in PAMs, and these mutations strongly decrease the targeting efficiency of bacterial spacers (Leenay et al., 2016; Sun et al., 2013).

In certain limits, average immunity can be related to the Morisita–Horn similarity index $\mathrm{\Psi }$Wolda, 1981 (see Appendix 1), where instead of comparing two populations from different timepoints or different regions, we compared bacteria and phage populations scaled by the immune benefit of CRISPR. By comparing our definition of average immunity to the Morisita–Horn index, we find that the Morisita overlap is constant across all parameters we study, which implies that the population diversity evolves to attain the highest possible overlap. Average immunity is also theoretically similar to the concept of distributed immunity introduced by Childs et al., 2014; Childs et al., 2012, the ‘immunity’ quantity calculated by Han et al., 2013, and the probability of an immune encounter used by Iranzo et al., 2013 to derive mean-field population predictions. The mathematical definition of average immunity is slightly different than these quantities, but the effect on populations is strikingly similar, and it is significant that several different models and different metrics find that a quantification of immunity is the important predictor for population-level outcomes.

We calculated average immunity in two experimental datasets and performed synthetic time-shift experiments by calculating average immunity between time-shifted populations of bacteria and phage. This effectively replicates explicit time-shift experiments (Richman et al., 2003; Frost et al., 2005; Moore et al., 2009; Blanquart and Gandon, 2013; Koskella, 2014; Chabas et al., 2016; Pyenson and Marraffini, 2020; Laanto et al., 2017; Common et al., 2019; Dewald-Wang et al., 2022). A similar time-shift calculation using sequence data was performed by Guillemet et al., 2021; they approached the question from the phage perspective and found that phages are most infectious against bacteria from the past and least infectious against bacteria from the future. Our results show that powerful insights into the dynamical immune state of a population can be obtained from sequencing data without explicitly performing time-shift experiments. We applied the same simple method to two very different experimental populations, suggesting that this analysis may be feasible in other time-series datasets, including from other natural microbial populations.

Our model predicts that an immune-evolving population in true steady state must eventually experience declining past immunity. CRISPR arrays are not infinite, and spacers must eventually be lost, even though this may happen on very long timescales. This pattern of declining past immunity is also in line with the expectations of fluctuating selection dynamics in which more rare genotypes are more fit (Gaba and Ebert, 2009; Hall et al., 2011; Blanquart and Gandon, 2013; Dewald-Wang et al., 2022). Previous theoretical work applied to vertebrate adaptive immunity also supports the existence of a peak in adaptive benefit in the recent past (Blanquart and Gandon, 2013; Nourmohammad et al., 2016), and some experimental time-shift studies have reported a peak in adaptation for bacteria-phage interactions (Hall et al., 2011; Koskella, 2014; Dewald-Wang et al., 2022 and for HIV-immune system interactions Blanquart and Gandon, 2013; Nourmohammad et al., 2016). We saw two qualitatively different time-shift curves in the experimental data we analysed: in the long-term laboratory coevolution data, bacteria were more immune to all past phages and less immune to all future phages (Figure 7C), while in the wastewater treatment plant data, bacteria were most immune to phages in their present context and less immune to phages in both the past and the future (Figure 7D). The time-shift pattern of the laboratory coevolution data was consistent with experimental time-shift data that is typically said to be indicative of arms-race dynamics (Richman et al., 2003; Frost et al., 2005; Moore et al., 2009; Blanquart and Gandon, 2013; Betts et al., 2018; Common et al., 2019), while the wastewater data may be consistent with a ‘zoomed-out’ view of the eventual decline in immunity in both the distant past and distant future (Nourmohammad et al., 2016; Guillemet et al., 2021). Indeed, the distinction between arms-race dynamics and fluctuating selection dynamics in time-shift experiments may be a question of the timescale being investigated (Gaba and Ebert, 2009; Blanquart and Gandon, 2013; Dewald-Wang et al., 2022). Our simulation results may qualitatively match both of these regimes depending on the timescale observed, which suggests that measuring very long-term coexistence is an area for further work in host-parasite coexistence experiments. In general, we expect that the timescale of memory length is related to the size of CRISPR arrays: if overall rates of spacer gain and loss are balanced, then we expect an individual spacer’s lifetime in the population to be proportional to array length. The factors that affect the timescale of memory length is an interesting area for further work.

One reason why past immunity did not unambiguously decline in our analysis of the laboratory coevolution data may be that spacer loss was not observed in the original experiment (Paez-Espino et al., 2015. This is also true of several other similar experiments with S. thermophilus Paez-Espino et al., 2013; Weissman et al., 2018), though spacer loss was directly observed coincident with spacer acquisition in Deveau et al., 2008. Spacer loss has also been observed in other experimental systems (Jiang et al., 2013; Rao et al., 2017; Deecker and Ensminger, 2020; Garrett, 2021), and indirect evidence of spacer loss through sequence comparisons has been reported for S. thermophilus (Horvath et al., 2008 and metagenomic data Held et al., 2010; Weinberger et al., 2012a). Spacers may also be functionally lost despite remaining in the genome, either from stochastic decoupling of RNA polymerase (Zoephel and Randau, 2013; Martynov et al., 2017; Soto-Perez et al., 2019 or from a dilution effect of competition for limited Cas protein complexes Martynov et al., 2017; Bradde et al., 2020; Garrett, 2021). The ‘loss’ rate in our model could also be taken to be gain and loss of the entire CRISPR system, which is known to happen (Delaney et al., 2012; Rollie et al., 2020), although this would likely happen on much slower timescales than the acquisition and loss we model. The qualitative shape and timescale of time-shift curves like the ones we presented may indirectly contain information about the timescale of spacer loss.

We found in our simulations that memory length decreased as spacer acquisition probability increased and as phage mutation rate decreased: if bacteria followed phages more closely (higher $\eta$) or if phages evolved more slowly (lower $\mu$), the peak in past immunity was closer to the present and memory was more quickly lost (Figure 7A, E, and F). A similar correspondence was derived in a theoretical model of optimal immune memory formation in vertebrates: more frequent sampling of the pathogen landscape (corresponding to higher spacer acquisition $\eta$ in our model) led to faster memory loss (Mayer et al., 2019). Notably, that result is the optimal update scheme for an immune system wishing to minimize the costs of infection; there is no explicit optimization in our model, yet we find a pattern consistent with a potentially optimal solution.

CRISPR adaptive immunity is mechanistically very different from the vertebrate adaptive immune system, yet there are several striking conceptual similarities between them and in the way past theoretical works have modelled vertebrate immunity. For example, phenomenological models of binding affinity between antigens and the immune system have used the degree of similarity between strings as a metric for interaction strength (Detours and Perelson, 1999; Detours and Perelson, 2000; Chao et al., 2005; Wang et al., 2015; Nourmohammad et al., 2016; Sachdeva et al., 2020), and more abstract models of the change in immunity after a virus mutation have also used string similarity to track the strength of immune coverage (Luksza and Lässig, 2014; Rouzine and Rozhnova, 2018; Yan et al., 2019). This is a simplification of a complex protein-protein interaction in the case of the vertebrate immune system (Altan-Bonnet et al., 2020), but it is nearly identical to the actual mechanism of CRISPR immunity, suggesting that lessons from models of CRISPR immunity may in turn be applicable to other forms of adaptive immunity.

Pathogen evolution has been extensively explored in models of vertebrate immunity, and here too we can draw conceptual parallels. We found that the speed of phage evolution depended sublinearly on phage population size in our model (Figure 6C), analogous to the weak positive dependence on population size described by Marchi et al., 2021. Adding cross-reactivity between spacer types in our model changed overall population outcomes as well as population dynamics. We observed a traveling wave regime in which new phage mutants ‘pull’ existing bacterial clones along, even if their matching clone is small or extinct (Figure 54). This travelling wave regime is qualitatively similar to that observed in theoretical models of ‘ballistic’ pathogen evolution in effective low-dimensional antigenic space (Yan et al., 2019; Marchi et al., 2019; Marchi et al., 2021). Similar clone population dynamics were observed in a model of B cell activation and memory generation in which existing memory B cells are reactivated in response to infection with a mutated virus that remains within its cross-reactivity radius (Chardès et al., 2022). We also observed splitting of lineages into divergent subtypes in our model, a phenomenon that has been studied in vertebrate pathogen evolution as well (Marchi et al., 2019; Marchi et al., 2021). We found that the average number of distinct clans was proportional to overall diversity but that increasing cross-reactivity decreased the average clan size (‘Number and size of clone clans’). These qualitative similarities highlight the usefulness of CRISPR adaptive immunity to understand population dynamics that also occur in the vertebrate immune system.

Selection in the vertebrate immune system happens at many distinct stages within individuals, including T cell receptor selection in the thymus (Camaglia et al., 2022 and affinity maturation of B cells following an immune challenge Chardès et al., 2022; Nourmohammad et al., 2019; Nourmohammad et al., 2016; Molari et al., 2020). Selection also happens at the population level as pathogens and individuals coevolve over time (Marchi et al., 2019; Yan et al., 2019, for instance, in the evolution of influenza Luksza and Lässig, 2014; Yan et al., 2019). In contrast, within-host dynamics are much simpler in bacteria, and the possibilities for within-host coevolution are especially limited in the case of virulent phages that kill their hosts after a single infection. Similarly, in the vertebrate immune system, extremely high immune diversity is possible in a single individual, whereas in CRISPR immunity, total diversity per individual is many orders of magnitude lower: most bacteria contain tens of spacers with a very few examples of a few hundred spacers (Pavlova et al., 2021). This reintroduces the question of how to think about immunity in microbial populations, whether as playing out on the level of individuals or on the level of the entire population. Finally, unlike the vertebrate adaptive immune system, CRISPR immunity is heritable, and so the optimal immune strategy for bacteria takes place on a timescale potentially much longer than an individual’s lifetime as well as on the level of the entire population Mayer et al., 2016. We have highlighted several analogies between CRISPR adaptive immunity and vertebrate adaptive immunity. Overall, the mechanistic simplicity and experimental tractability of CRISPR immunity mean that abstract ideas from vertebrate immunity can be linked to measurable features of CRISPR immunity, and insights from CRISPR immunity may be useful in understanding vertebrate immunity as well.

In some CRISPR systems, an imprecise match between a spacer and protospacer results in a much higher spacer acquisition rate, a phenomenon called priming (Rao et al., 2017). Including cross-reactivity in our model is a way to explore the effect of primed acquisition. However, we modelled cross-reactivity as a change in phage infection success probability, whereas in primed acquisition it is the spacer acquisition rate that changes. Because bacteria must lose their spacer before gaining a new spacer in our model, including primed acquisition directly would require expanding bacterial arrays to at least two spacers. The acquisition probability $\eta$ could then become a function of protospacer-spacer similarity to directly model primed acquisition.

The mechanism of selection for immune variants in CRISPR adaptive immunity is a matter of debate. We observed several features that suggest the presence of negative frequency-dependent selection in our model: (a) decreased clone growth rates as clone size increased, (b) cyclic clone size dynamics as immune pairs oscillate in size, and (c) continual turnover in spacer types over time. This is consistent with recent experimental observations of negative frequency-dependent selection in host-pathogen systems (Betts et al., 2018; Guillemet et al., 2021). Our model also contains the base assumption that phages that can successfully infect dominant bacterial strains will experience positive selection, generating ‘kill-the-winner‘ dynamics that lead to negative frequency-dependent selection for bacteria (Weinberger et al., 2012a). Our model does not capture true host-pathogen coevolution because bacteria are not able to evolve beyond ‘following’ the mutational landscape of phages by acquiring spacers. In true coevolution, bacterial genomic mutations, apart from the ability to acquire spacers, are also drivers of evolution and sources of selection pressure on phages.

We observed large population size fluctuations when average immunity reached a particular value of approximately 10% for the parameters we used (Figure 31). This critical value is determined by a balance of the rate of spacer acquisition and spacer loss that leads to exactly equal fitness for both bacteria with spacers and bacteria without spacers (see ‘Dominant balance approximations for $\nu$’). Despite the presence of spacers with non-zero effectiveness, near this critical point the population appears to behave as if there is no CRISPR immunity at all. This has two interesting results: first, the probability of phage establishment drops to zero for some simulations near this critical point (Figure 3—figure supplement 2), since if CRISPR has no influence then mutant phages do not experience strong positive selection. Second, population sizes become unstable near the transition: outcomes with quite different total bacteria population size and fraction of bacteria with spacers are possible (Figures 30 and 31). This critical point will be interesting to explore further as a possible phase transition.

We found that diversity in our model depended on spacer acquisition rate, spacer effectiveness, and spacer mutation rate. These parameters are all known to vary in natural systems, and some may also be experimentally manipulated. CRISPR systems display a range of spacer acquisition rates even within the same system (Heler et al., 2019), and both spacer acquisition rate and spacer effectiveness may be controlled by the level of expression of Cas proteins, which has been found to be connected to the quorum sensing pathway (Høyland-Kroghsbo et al., 2017). Additionally, recent work has found that spacer acquisition is more efficient at slow bacterial growth rate because phages also grow more slowly, giving bacteria more time to acquire spacers (Høyland-Kroghsbo et al., 2018; Dimitriu et al., 2022). This change in bacterial growth rate can be controlled by nutrient concentration (Dimitriu et al., 2022; Payne et al., 2018, temperature Høyland-Kroghsbo et al., 2018, or by bacteriostatic antibiotics which slow bacterial growth without killing Dimitriu et al., 2022). Phage mutations rates also vary and may be experimentally modified (Kysela and Turner, 2007). These findings all represent possible experimental manipulations of parameters relevant for CRISPR immunity; for example, given these findings we predict that a bacteria-phage population exposed to bacteriostatic antibiotics would evolve higher CRISPR spacer diversity than one without antibiotics.

In this section, we describe the details of our phenomenological model of bacteria and phages, providing more detail for Figure 1A. Sections ‘Reactions’ and ‘Master equation’ list the stochastic interactions we simulate, and ‘Clone size master equation’ lists master equations that are used to derive mean-field equations for this dynamical system (Total population size) and clone size distributions for bacteria and phages (Appendix 2).

We model bacteria and phages interacting in a chemostat. The populations we track are nutrient concentration ($C$), phages ($n_V$), and bacteria ($n_B$). Bacteria can either have no spacer ($n_B^0$) or a spacer of type $i$ ($n_B^i$), and phages can have a single protospacer of type $j$ ($n_V^j$). Nutrients flow in at concentration C₀ with rate $F$, and all species flow out with rate $F$. The total number of bacteria with a spacer is $n_B^s$ and the total number of bacteria is $n_B$. Bacteria grow at rate $gC$. With rate $\alpha$, a phage interacts with a bacterium. With probability $p_V$, the phage will kill bacteria without a matching spacer and produce a burst of new phages with size $B$, while for bacteria with a matching spacer that probability is reduced to $p_V^s=(1-e)p_V$ ($0\le e\le 1$). Bacteria without spacers that survive an attack have a chance to acquire a spacer with probability $\eta$. Bacteria with a spacer lose their spacer at rate $r$. Phages are limited to a single protospacer, and phages can mutate to a new protospacer type with probability $\mu L$, where $\mu$ is the per-base mutation rate and $L$ is the protospacer length ($L=30$ in all simulations). Parameter descriptions and default values are shown in Table 3. See refs Weissman et al., 2021; Gurney et al., 2019 for recent examples of other chemostat-based models.

This section describes how simulations were performed, shows basic simulation results, describes parameter choices (‘Parameter values’), and explores the parameter dependence of stochastic population extinction (‘Stochastic population extinction’) that relates to results in ‘Bacteria and phages dynamically coexist and coevolve’.

Simulations were written in Python and performed on the Béluga and Niagara supercomputers at the SciNet HPC Consortium (Ponce et al., 2019) and on a local server in our group. Our simulation code can be found at https://github.com/mbonsma/CRISPR-dynamics-model (copy archived at Bonsma-Fisher and Goyal, 2023). We used the tau-leaping method, a fast approximation for Gillespie simulation Cao et al., 2006 and compared with full Gillespie simulations for some cases. Gillespie and tau-leaping methods showed similar dynamics and good agreement for total population sizes (Figure 11), total spacer diversity (Figure 12), mean phage and bacteria clone sizes (Figure 13), and produced the same qualitative behaviour for individual spacer types (Figure 14).

Figure 11

Download asset Open asset

Figure 12

Download asset Open asset

Figure 13

Download asset Open asset

Figure 14

Download asset Open asset

We initialized five simulations for each parameter combination for a total of approximately 19,800 simulations, not including simulations with cross-reactivity. A small subset of simulation parameters were run six times instead of five. Not all simulations were successfully completed either due to errors while running or because their running time was exceedingly long. Of the initialized simulations, approximately 10,000 were completed and included in analysis. Simulations with very low $\mu$ tended to either have no phage establishments (at low C₀) or very long running times (at high C₀), and simulations with very high $\mu$ tended to have higher diversity and longer running times in general. Unless otherwise noted, plots with simulation averages are an average across three to six simulations with the same parameters. We set simulations to run for a fixed number of bacterial generations intended to be long enough to allow the system to reach steady state and remain there for a long time to generate statistics (‘Steady state’).

We model spacers and protospacers as a binary sequence of length $L$ as in Han et al., 2013; Han and Deem, 2017. When a phage reproduces and creates a burst of $B$ phages, we draw $BL$ numbers from a binomial distribution with probability $\mu$. Successes in this draw designate bits that will mutate (flip) in the newly created phages. It is possible but very rare for more than one mutation to happen in the same new phage: this occurs with probability $1-{\text{e}}^{-\mu L}-\mu L{\text{e}}^{-\mu L}\approx {(\mu L)}^2$ whereas a single mutation occurs with probability $\mu L{\text{e}}^{-\mu L}\approx \mu L$. Multiple mutations occur approximately $\mu L$ times as often as single mutations; for the largest value of $\mu$ we use, they represent 0.3% of events.

Our simulation code can be found on GitHub.