Figures and data in Dynamics of immune memory and learning in bacterial communities

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: April 20, 2023 (This version)
Accepted Manuscript published: January 16, 2023 (Go to version)
Accepted: January 15, 2023
Received: July 7, 2022
Preprint posted: July 7, 2022 (Go to version)

Figures
Tables
Additional files

102 figures, 4 tables and 1 additional file

Figures

Figure 1 with 7 supplements

Download asset Open asset

Model description.

(A) We model bacteria and phages interacting in a well-mixed vessel. We track nutrient concentration, phage population size ( $n_{V}$ ), and bacteria population size ( $n_{B}$ ). Bacteria can either have no spacer ( $n_{B}^{0}$ ) or a spacer of type $i$ ( $n_{B}^{i}$ , $\sum_{i} n_{B}^{i} = n_{B}^{s}$ ), and phages can have a single protospacer of type $j$ ( $n_{V}^{j}$ ). With rate $α$ , a phage interacts with a bacterium. If the bacterium does not have a matching spacer, the phage kills with probability $p_{V}$ and produces a burst of $B$ new phages, while for bacteria with a matching spacer that probability is reduced to $p_{V}^{s} = p_{V} (1 - e)$ , $0 \leq e \leq 1$ . Bacteria without spacers that survive an attack have a chance to acquire a spacer with probability $η$ , and bacteria with spacers lose them at rate $r$ . Lower inset: average immunity is the weighted average pairwise immunity between spacer-containing bacteria and phages, given by $1 - \frac{\sum_{i, j} n_{B}^{i} n_{V}^{j} p_{V} (i, j)}{p_{V} \sum_{i, j} n_{B}^{i} n_{V}^{j}}$ . The probability of a phage with protospacer $j$ successfully infecting a bacterium with spacer $i$ is $p_{V} (i, j)$ . (B) Three time points in a typical simulation with $C_{0} = 10^{4}$ , $e = 0.95$ , $η = 10^{- 4}$ , and $μ = 10^{- 5}$ . Coloured circles represent unique protospacer or spacer sequences; shared sequences are shown with the same colour. The size of each circle is proportional to clone size, and new mutants are shown radially more distant from the centre. (C) Ten individual clone trajectories vs simulation time for phages (top) and bacteria (bottom). The mean clone size is shown with a horizontal dashed line. (D) Total phage, bacteria, and nutrient concentration as a function of phage success probability $p_{V}$ . Markers show an average over five independent simulations for different values of $p_{V}$ with $C_{0} = 10^{4}, η = 10^{- 3}, e = 0.95$ , and $μ = 10^{- 7}$ . Solid lines show theoretical predictions for different constant values of effective $e$ . As $p_{V}$ decreases, phages go extinct at a critical value given by $A = 1$ , where $A = \frac{(B p_{V} - 1) (1 - f) α}{f g}$ . (E) Total phage and bacteria population size as a function of average bacterial immunity to phages. Colours indicate the fraction of simulations in which phage or bacteria go extinct before a set endpoint. Solid lines show the mean-field prediction. Error bars are the standard deviation across three or more independent simulations.

Figure 1—figure supplement 1

Download asset Open asset

Probability of stochastic extinction at low spacer acquisition.

Probability of extinction in four or more simulations with the same parameters vs. mean phage population size (left) and mean bacteria population size (right) for the lowest value of spacer acquisition probability ( $η = 10^{- 5}$ ). Colours indicate mutation rate $μ$ and shapes indicate CRISPR effectiveness $e$ .

Figure 1—figure supplement 2

Download asset Open asset

Probability of stochastic extinction at high spacer acquisition.

Probability of extinction in four or more simulations with the same parameters vs. mean phage population size (left) and mean bacteria population size (right) for the highest value of spacer acquisition probability ( $η = 10^{- 2}$ ). Colours indicate mutation rate $μ$ and shapes indicate CRISPR effectiveness $e$ . Error bars are the standard deviation across three or more independent simulations.

Figure 1—figure supplement 3

Download asset Open asset

Time to extinction for phages vs. mutation rate.

Mean time for the phage population to go extinct vs. phage mutation rate $μ$ across simulations where at least one simulation experienced phage extinction. The darkness of each point indicates the fraction of simulations that went extinct with darkest colours representing all simulations extinct. Simulations are shown for $η = 10^{- 2}$ . Error bars are the standard deviation across three or more independent simulations.

Figure 1—figure supplement 4

Download asset Open asset

Time to extinction for bacteria vs. mutation rate.

Mean time for the bacteria population to go extinct vs. phage mutation rate $μ$ across simulations where at least one simulation experienced bacteria extinction. The darkness of each point indicates the fraction of simulations that went extinct with darkest colours representing all simulations extinct. Error bars are the standard deviation across three or more independent simulations.

Figure 1—figure supplement 5

Download asset Open asset

Time to extinction for phages for different initial diversity and low spacer acquisition.

Mean time for the phage population to go extinct vs. phage mutation rate $μ$ across simulations where at least one simulation experienced phage extinction. The darkness of each point indicates the fraction of simulations that went extinct with darkest colours representing all simulations extinct. Colours represent number of initial phage clones $m_{i n i t}$ . Simulations are shown for $η = 10^{- 5}$ and $C_{0} = 1000$ . Error bars are the standard deviation across three or more independent simulations.

Figure 1—figure supplement 6

Download asset Open asset

Time to extinction for phages for different initial diversity and high spacer acquisition.

Mean time for the phage population to go extinct vs. phage mutation rate $μ$ across simulations where at least one simulation experienced phage extinction. The darkness of each point indicates the fraction of simulations that went extinct with darkest colours representing all simulations extinct. Colours represent number of initial phage clones $m_{i n i t}$ . Simulations are shown for $η = 10^{- 2}$ and $C_{0} = 1000$ . Error bars are the standard deviation across three or more independent simulations.

Figure 1—video 1

Download asset

posterframe for video — An animation of a typical simulation of bacteria and phages interacting with CRISPR immunity.

Each bacteria can have a single spacer, indicated by colour, and each phage has a single protospacer sequence. Matching sequences are shown with the same colour. The size of each circle is proportional to clone size, and new mutants are shown radially more distant from the centre. New phages appear by mutation over time; many of them go extinct, but some establish and grow large, and then bacteria acquire matching spacers.

Figure 2 with 4 supplements

Download asset Open asset

Diversity depends sub-linearly on parameters.

(A) Bacteria and phage clone size distributions normalized to the measured mean clone size for $C_{0} = 10^{5}$ , $μ = 3 \times 10^{- 7}$ , and $e = 0.95$ . As $η$ increases, both clone size distributions become more sharply peaked. (B) The mean number of bacterial clones depends only on a combined parameter in the limit of small average immunity (generally coinciding with high C₀). (Inset) The mean number of bacterial clones can be predicted by numerically solving Equation 1 for $m$ . The two lowest values of $η$ are shown with lighter shading. Error bars are the standard deviation across three or more independent simulations.

Figure 2—figure supplement 1

Download asset Open asset

Predicted diversity grouped by spacer effectiveness.

Simulation mean number of bacteria clones ( $m$ ) vs. theoretical prediction for $m$ given by numerically solving Equation 1, broken down by different values of $e$ . Error bars are the standard deviation across three or more independent simulations.

Figure 2—figure supplement 2

Download asset Open asset

Predicted diversity grouped by mutation rate.

Simulation mean number of bacteria clones ( $m$ ) vs. theoretical prediction for $m$ given by numerically solving Equation 1, broken down by different values of $μ$ . Error bars are the standard deviation across three or more independent simulations.

Figure 2—figure supplement 3

Download asset Open asset

Clone size histograms by spacer effectiveness.

Bacteria and phage clone size distributions normalized to the measured mean clone size for $C_{0} = 10^{5}$ , $μ = 3 \times 10^{- 7}$ , and $η = 0.001$ .

Figure 2—figure supplement 4

Download asset Open asset

Clone size histograms by phage mutation rate.

Bacteria and phage clone size distributions normalized to the measured mean clone size for $C_{0} = 10^{5}$ , $e = 0.95$ , and $η = 0.001$ .

Figure 3 with 3 supplements

Download asset Open asset

The fate of individual clones.

(A) Phage and bacteria coevolve in two timescale-separated regimes characterized by phage clone fitness. Average phage and bacteria clone size vs. time since phage mutation (right axis), and average clone growth rate vs. time since phage mutation (left axis). Markers show the average over all clone trajectories after steady state from six simulations with the same parameters. (B) Histograms of individual clone fitness grouped by time since phage mutation. Phage clones initially have fitness gt₀, but rapidly most clones reach neutral growth (fitness $\approx 0$ ). Bacteria clones also follow suit, initially having fitness gt₀ and rapidly reaching 0 fitness on average. Because spacer acquisition for a clone only happens after that clone is created by phage mutation, the top-right panel of (B) is empty at the earliest time point following phage mutation. Individual clone trajectories are highly variable. (C) Probability of phage clone establishment vs. average immunity. Clones are considered established in simulations when they reach the mean clone size. Equation 3 with $ν = 1$ is shown in green and with $ν$ given by Equation 4 in blue. In (**A, B**), $C_{0} = 10^{4}$ , $e = 0.95$ , $η = 10^{- 3}$ , and $μ = 3 \times 10^{- 6}$ . Error bars are the standard deviation across three or more independent simulations.

Figure 3—figure supplement 1

Download asset Open asset

Average clone sizes and fitness over time in a simulation.

Average phage and bacteria clone size vs. time (right vertical axis, purple and green markers) and average phage clone growth rate vs. time (left vertical axis, orange markers). Markers are the average over all clone trajectories after steady state in a single simulation. Phage clones appear at size 1 and grow until they reach the deterministic mean phage clone size on average. Once a phage clone becomes large, bacteria encounter it often enough to acquire a matching spacer, and bacteria clones grow until they reach their deterministic mean clone size. New phage mutants have a selective advantage (fitness gt₀, positive growth rate) until they reach the deterministic mean clone size, at which point they evolve neutrally (fitness $\approx 0$ ). Individual clone trajectories are highly variable, leading to a large standard deviation on the mean phage clone fitness and a fitness trend which is only evident on average. The predicted clone size is piecewise-defined as the theoretical clone trajectory until the theoretical trajectory reaches the predicted mean clone size.

Figure 3—figure supplement 2

Download asset Open asset

Approximations for phage establishment probability.

Phage establishment probability vs. effective $e$ . Markers are simulation results; error bars are standard deviation across three or more independent simulations. Black lines are the numerical theoretical solution. Blue lines are the approximate solution given by Equation 54 (with the full theoretical predicted value of $m$ ), and green lines are the same but with $e = 0$ . All values of $e$ collapse on the same $η$ lines; the establishment probability depends only on $e / m$ and not on $e$ by itself. Interestingly, the $e = 0$ approximation is close to the small $e$ approximation at high $η$ . This implies that high $e$ matters more at small $η$ . For some simulations near the transition from low to high $η$ , the probability of establishment is zero, which is why the grey connecting lines drop below the x-axis near effective $e \approx 10^{- 1}$ .

Figure 3—figure supplement 3

Download asset Open asset

Theoretical phage establishment probability with approximations.

Theoretical phage establishment probability vs. $e / m$ . Solid lines are the numerical theoretical solution. Light dashed lines (‘small $e_{e f f}$ ’) are the approximate solution given by Equation 54 (with the full theoretical predicted value of $m$ ), heavy dashed lines are the same but with $e = 0$ . All values of $e$ collapse on the same $η$ lines; the establishment probability depends only on $e / m$ and not on $e$ by itself. The inflection point in $P_{e s t}$ (dashed black line) is given by Equation 58.

Figure 4 with 4 supplements

Download asset Open asset

Cross-reactivity leads to ‘spindly’ phylogenies and regime switching.

Phage clone phylogenies for four simulations with different cross-reactivities: no cross-reactivity (A), exponential cross-reactivity with $θ = 4$ (B), and step-function cross-reactivity with $θ = 1$ (C) and $θ = 2$ (D). All simulations share all other parameters: $C_{0} = 10^{4}, η = 10^{- 4}, μ = 10^{- 6}, e = 0.95$ . Phage clones are plotted at the first time they pass a population size of 2 to remove clutter from many new mutations destined for extinction, and the size of each circle is logarithmically proportional to the maximum size reached by that clone. Colours indicate the time of extinction of each clone. For each simulation with cross-reactivity, the left inset shows phage (top) and bacteria (bottom) clone sizes over time; colours indicate unique clone identities. Coloured rectangles above insets in (C) and (D) correspond to the dominant clone at each time. Dominant clone identities are offset by $θ$ (vertical dashed line for visual aid).

Figure 4—video 1

Download asset

Figure 4—video 2

Download asset

Figure 4—video 3

Download asset

Figure 4—video 4

Download asset

Figure 5 with 1 supplement

Download asset Open asset

Average immunity underlies population outcomes.

(A) Probability of phage clone establishment vs. average immunity for different amounts and types of cross-reactivity. No cross-reactivity ( $θ = 0$ ) is shown as black stars, exponential cross-reactivity in red, and step-function cross-reactivity in blue. Simulation averages are shown for $η = 10^{- 4}$ and $μ = 10^{- 6}$ . Error bars are the standard deviation across three or more independent simulations and are shown in both x directions and the positive y direction. (**B, C**) Total phage (purple) and total bacteria (teal) average population sizes vs. the mean number of bacterial clones $m$ (B) and vs. average bacterial immunity (C) for $η = 10^{- 4}$ . Each point is an average at steady state over three or more independent simulations with the same parameters; error bars are standard deviation. Total sizes are scaled by the initial nutrient concentration C₀. Lighter colours indicate stronger cross-reactivity, marker shapes match legends in (A) and (B). Solid lines are the predicted total population size given by solving Equations 13–17 and using the approximation effective $e \approx e / m$ in (B) and the measured average immunity for effective $e$ in (C).

Figure 5—figure supplement 1

Download asset Open asset

Total population sizes predicted by average immunity in a simulation with cross-reactivity.

Total phage (top) and total bacteria (bottom) in a simulation with cross-reactivity (step function CRISPR effectiveness with $θ = 1$ ). The dashed black line uses the measured value of average immunity from the simulation at each time point to predict population sizes using the solutions to the system of Equations 13–17. This simulation has one initial phage clone and parameters $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 10^{- 4}$ , and $e = 0.95$ .

Figure 6 with 5 supplements

Download asset Open asset

Phage evolution and spacer turnover.

(A) Principal Component Analysis (PCA) decomposition of phage and bacteria clone abundances for a simulation with $C_{0} = 10^{4}$ , $e = 0.95$ , $η = 10^{- 4}$ , and $μ = 10^{- 5}$ . Clone abundances are normalized at each time point, then PCA is performed for the entire phage time series over $\approx 4000$ generations (four times the mean extinction time for phage clones). Bacteria and phage clone abundances are transformed into the PCA coordinates; colours indicate simulation time. Five time points are highlighted in progressively lighter shades of red for emphasis. (B) Phage genomic speed of evolution vs. mean large phage clone time to extinction. The phage speed is the weighted average genomic distance between the phage population at the end of the simulation and the phage population at an earlier time, divided by the time interval. The dashed line is $y = \frac{1}{x}$ . (C) The speed of evolution increases as spacer effectiveness $e$ , spacer acquisition probability $η$ , and phage mutation rate $μ$ increase. The dashed line shows an approximate theoretical calculation (assuming speed = 1/time to extinction) which captures the trend across a wide range of parameters. Error bars in (B) and (C) are the standard deviation across three or more independent simulations and are shown in the positive direction only. (D) Spacer turnover as a function of time delay for four simulations with $C_{0} = 10^{4}$ , $e = 0.95$ , and $μ = 10^{- 5}$ . The fraction of bacterial clones remaining is the fraction of clones that were present at time $t$ that are still present at time $t +$ delay. Solid lines are an average across steady-state for each value of the time delay; shaded regions are the standard deviation. (**E–G**) Spacer-type turnover calculated as in (D) using experimental data from Paez-Espino et al., 2015 (E), metagenomic data sampled from groundwater from Burstein et al., 2016 (F), and metagenomic data sampled from a wastewater treatment plant from Guerrero et al., 2021a (G). Experimental time points are interpolated to the minimum sampling interval to allow averaging across the experiment.

Figure 6—figure supplement 1

Download asset Open asset

Clone size PCA for a simulation with exponential cross-reactivity.

PCA decomposition of phage and bacteria clone abundances for a simulation with exponential cross-reactivity and $C_{0} = 10^{4}$ , $e = 0.95$ ,, $η = 10^{- 4}$ and. $μ = 10^{- 6}$ Clone abundances are normalized at each time point, then PCA is performed for the entire phage time series over $\approx 500$ generations (four times the mean extinction time for phage clones). Bacteria and phage clone abundances are transformed into the PCA coordinates; colours indicate simulation time. Five time points are highlighted in progressively lighter shades of red for emphasis.

Figure 6—figure supplement 2

Download asset Open asset

Genetic distance vs. diversity and time to extinction.

Phage mutational distance reached in simulations divided by the number of phage establishments vs. diversity (left) and phage mutational distance per generation vs. mean phage time to extinction (right). Error bars are the standard deviation across three or more independent simulations and are shown in the positive direction only.

Figure 6—figure supplement 3

Download asset Open asset

Spacer turnover by spacer effectiveness.

Spacer turnover as a function of time delay for four simulations with $C_{0} = 10^{4}$ , $e t a = 0.001$ , and $μ = 10^{- 5}$ . The fraction of bacterial clones remaining is the fraction of clones that were present at time $t$ that are still present at time $t +$ delay. Solid lines are an average across steady state for each value of the time delay; shaded regions are the standard deviation.

Figure 6—figure supplement 4

Download asset Open asset

Spacer turnover by phage mutation rate.

Spacer turnover as a function of time delay for four simulations with $C_{0} = 10^{4}$ , $e = 0.95$ , and $e t a = 0.001$ . The fraction of bacterial clones remaining is the fraction of clones that were present at time $t$ that are still present at time $t +$ delay. Solid lines are an average across steady state for each value of the time delay; shaded regions are the standard deviation.

Figure 6—video 1

Download asset

Figure 7 with 12 supplements

Download asset Open asset

Quantifying immune memory in data.

(**A, B**) Average immunity of bacteria against phage for four simulations with different values of $η$ as a function of time shift. Solid lines are an average across steady state for each value of the time shift; shaded regions are the standard deviation. Average immunity peaks in the recent past (A, indicated by $τ^{*}$ ) with a negative slope through zero delay (A, black dashed line) and decays to zero at long delays in the past or future (B). For all simulations $C_{0} = 10^{4}$ , $μ = 10^{- 5}$ , and $e = 0.95$ . (**C, D**) Average overlap between bacterial spacer and phage protospacer types using data from a lab experiment with *S. thermophilus* and phage from Paez-Espino et al., 2015 (C) and data from a wastewater treatment plant sampled over 3 years from Guerrero et al., 2021a (D). Spacer types are grouped by 85% similarity, and shaded region is standard deviation across averaged data. Base average immunity values were multiplied by the average number of protospacers corresponding to the *S. thermophilus* CRISPR system (C) and the *Gordonia* CRISPR systems (D) to account for multiple potential protospacer targets per phage. In (D), we compared two time shifts with zero delay average immunity using a Wilcoxon signed-rank test: $p = 0.27$ for lower past immunity at 500 days, $p = 0.008$ for lower past immunity at 200 days, $p = 0.001$ for lower future immunity at 500 days, and $p = 0.003$ for lower future immunity at 200 days. (**E, F**) The position of the peak in past immunity for simulated data vs. spacer acquisition probability $η$ (E) and phage mutation rate $μ$ (F). The peak position is the time shift value for which the curves in (A) are largest, indicated by $τ^{*}$ . Error bars are the standard deviation across three or more independent simulations.

Figure 7—figure supplement 1

Download asset Open asset

Time-shifted average immunity for experimental coevolution data.

Time shifted overlap for experimental coevolution data with the first time point and last two time points removed. Interpolation spacing is 3 days, the smallest interval between remaining time points. Time intervals in days were multiplied by 6.64, the estimated number of bacterial generations per day assuming exponential growth with 100:1 serial dilutions. Protospacers are included if they have a perfect PAM, and all wild-type spacers are included.

Figure 7—figure supplement 2

Download asset Open asset

Time-shifted average immunity for experimental coevolution data with partial PAM matches included.

Time-shifted overlap for experimental coevolution data with the first time point and last two time points removed. Interpolation spacing is 3 days, the smallest interval between remaining time points. Time intervals in days were multiplied by 6.64, the estimated number of bacterial generations per day assuming exponential growth with 100:1 serial dilutions. Protospacers are included if they have a perfect PAM or partially perfect PAM, and all wild-type spacers are included.

Figure 7—figure supplement 3

Download asset Open asset

Time-shifted average immunity for experimental coevolution data with all sequences regardless of PAM.

Time-shifted overlap for experimental coevolution data with the first time point and last two time points removed. Interpolation spacing is 3 days, the smallest interval between remaining time points. Time intervals in days were multiplied by 6.64, the estimated number of bacterial generations per day assuming exponential growth with 100:1 serial dilutions. All potential protospacers are included regardless of PAM sequence, and all wild-type spacers are included.

Figure 7—figure supplement 4

Download asset Open asset

Time-shifted average immunity for experimental coevolution data with data trimmed from start and end.

Time-shifted overlap for experimental coevolution data with increasing numbers of time points from the start and end removed. Interpolation spacing is the smallest interval between remaining time points. Only perfect-PAM protospacers are included; all wild-type spacers are included.

Figure 7—figure supplement 5

Download asset Open asset

Time-shifted average immunity for wastewater with data trimmed from start and end.

Time-shifted average immunity for spacers and protospacers from wastewater grouped at four different similarity thresholds. Points are trimmed from the start and end of the time series in each panel as described in the overlay text. Raw values are multiplied by 1956, the average number of protospacers with the GTT PAM from the phage DC-56 and DS-92 genomes. Only perfect-PAM protospacers are included; all wild-type spacers are included.

Figure 7—figure supplement 6

Download asset Open asset

Bootstrapped average immunity with all points randomly shuffled.

Bootstrapped control: time-shifted average immunity for wastewater after randomly shuffling the interpolated clone abundances separately for bacteria and phage.

Figure 7—figure supplement 7

Download asset Open asset

Bootstrapped average immunity with pairs of bacteria and phage clone sizes shuffled.

Bootstrapped control: time-shifted average immunity for wastewater after randomly shuffling pairs of interpolated bacteria and phage clone sizes such that time point matching is maintained between bacteria and phages.

Figure 7—figure supplement 8

Download asset Open asset

Time-shifted average immunity for wastewater with data trimmed from the end.

Time-shifted average immunity for spacers and protospacers grouped at an 85% similarity threshold. 16 time points are removed from the end of the data series in order to remove the region with zero protospacer counts. Raw values are multiplied by 1956, the average number of protospacers with the GTT PAM from the phage DC-56 and DS-92 genomes. Only perfect-PAM protospacers are included; all wild-type spacers are included.

Figure 7—figure supplement 9

Download asset Open asset

Time-shifted average immunity for wastewater with data trimmed from the start.

Time-shifted average immunity for spacers and protospacers grouped at an 85% similarity threshold. 21 time points are removed from the beginning of the data series to remove the region with anomalously high average immunity. Raw values are multiplied by 1956, the average number of protospacers with the GTT PAM from the phage DC-56 and DS-92 genomes. Only perfect-PAM protospacers are included; all wild-type spacers are included.

Figure 7—figure supplement 10

Download asset Open asset

Time-shifted average immunity for wastewater with data trimmed from the start and end.

Time-shifted average immunity for spacers and protospacers grouped at an 85% similarity threshold. 21 time points are removed from the beginning of the data series to remove the region with anomalously high average immunity, and 16 time points are removed from the end of the data series in order to remove the region with zero protospacer counts. Raw values are multiplied by 1956, the average number of protospacers with the GTT PAM from the phage DC-56 and DS-92 genomes. Only perfect-PAM protospacers are included; all wild-type spacers are included.

Figure 7—figure supplement 11

Download asset Open asset

Time-shifted average immunity by spacer effectiveness.

Average immunity of bacteria against phage for four simulations with different values of $e$ as a function of time shift. Solid lines are an average across steady-state for each value of the time shift; shaded regions are the standard deviation. A short-time view is shown on the left and a long-time view on the right. For all simulations $C_{0} = 10^{4}$ , $e t a = 0.001$ , and $μ = 10^{- 5}$ .

Figure 7—figure supplement 12

Download asset Open asset

Time-shifted average immunity by phage mutation rate.

Average immunity of bacteria against phage for four simulations with different values of $μ$ as a function of time shift. Solid lines are an average across steady state for each value of the time shift; shaded regions are the standard deviation. A short-time view is shown on the left and a long-time view on the right. For all simulations $C_{0} = 10^{4}$ , $e = 0.95$ , and $e t a = 0.001$ .

Figure 8

Download asset Open asset

Total phage, total bacteria (left), and mean number of bacterial clones (right) vs. simulation time for five simulations with $μ = 10^{- 6}$ , $e = 0.5$ , $η = 0.001$ , and C₀ ranging from 300 (top row) to 30,000 (bottom row).

Total population sizes equilibrate very quickly, but the total number of clones can take longer at large population sizes (high C₀). The time constants inset are a measure of how quickly we expect each mean-field quantity to equilibrate: $n_{V}$ time constant is the inverse growth rate of the total phage population ( $1 / (- F - α n_{B} (B p_{V} - 1) - α B p_{V} n_{B}^{s} e / m)$ ) and the extinction time constant is the mean time to extinction for large phage clones (Equation 171), a measure of the rate of turnover of the number of clones.

Figure 9

Download asset Open asset

Mean number of bacterial clones after $t = t_{s s}$ bacterial generations vs. initial number of phage clones for $e = 0.8$ , $η = 10^{- 3}$ .

The mean is an average of 15 evenly spaced points from $t = t_{s s}$ to $t = 5 t_{s s}$ bacterial generations. Error bars are the standard deviation across three or more independent simulations.

Figure 10

Download asset Open asset

Phage (A) and bacteria (B) mean clone size in a simulation, either conditioned on survival (blue circles) or including extinct clones (orange circles).

Theoretical predictions are plotted as solid lines: the time-dependent numerical solution to Equations 21 and 22 in green, the same solution divided by the phage clone probability of survival in red, and a one-dimensional solution to Equation 21 in (A) and 22 in (B). Equations 25 and 26 are black dashed lines. An alternate simulation mean clone size is plotted for bacteria (brown circles) in which each clone trajectory is stacked based on the bacterial acquisition time and averaged across trajectories, conditioned on survival. Simulation parameters are $C_{0} = 10000.0$ , $μ = 10^{- 5}$ , $η = 0.001$ , and $e = 0.95$ .

Figure 11

Download asset Open asset

Total phage ( $n_{V}$ ) and total bacteria ( $n_{B}$ ) as a function of time for a Gillespie simulation and a tau-leaping simulation with the same parameters.

Total phage is shown at early times (A) and late times (B), and total bacteria at early times (C) and late times (D). The simulation parameters are $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 0.001$ , and $e = 0.95$ . Early time dynamics differ slightly in this example, but the long-time behaviour and steady-state values are similar.

Figure 12

Download asset Open asset

Number of bacterial clones ( $m$ ) vs. simulation time for three sets of simulations, each beginning with 1, 10, or 50 phage clones.

Gillespie simulations are dashed blue lines, and tau-leaping simulations are solid orange and red lines. The simulation parameters are $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 0.001$ , and $e = 0.95$ .

Figure 13

Download asset Open asset

Average population size of phage clones (solid lines) and bacterial clones (dashed lines) in a Gillespie and tau-leaping simulation with the same parameters.

The simulation parameters are $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 0.001$ , and $e = 0.95$ .

Figure 14

Download asset Open asset

10 spacer trajectories for a Gillespie simulation (A) and tau-leaping simulation (B).

The first 10 trajectories that surpass $n_{V}^{i} = 1000$ are shown. The simulation parameters are $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 0.001$ , and $e = 0.95$ .

Figure 15

Download asset Open asset

Phage clone size distribution from 15 combined time points for a simulation with the parameters $C_{0} = 10^{4}$ , $e = 0.95$ , $η = 0.01$ , and $μ = 3 \times 10^{- 6}$ .

The blue points are the values of the full normalized phage clone size histogram with a bin width of 1500. The orange line is given by $P_{n}^{large}$ in Equation 28 smoothed with a running average of window size 3000. Both distributions are scaled by the total number of phage clones.

Figure 16

Download asset Open asset

Mean phage clone size at the time of first spacer acquisition vs. the deterministic mean phage clone size.

For each phage clone trajectory, the clone size at the time of first spacer acquisition is recorded and these are averaged across each simulation. Error bars are the standard deviation across three or more independent simulations.

Figure 17

Download asset Open asset

Mean time to extinction for phage clones vs. the timescale of bacteria spacer acquisition given by $1 / D$ where $D = α η (1 - p_{V}) n_{V}^{i}^{*} n_{B}^{0} g C_{0}$ .

Points outlined in red are simulations where the ratio of large phage clones to bacterial clones exceeds 1.2. Phage clones experience clonal interference at low $η$ : they go extinct faster than bacteria acquire spacers.

Figure 18

Download asset Open asset

Ratio of average phage clone initial fitness to the average phage clone fitness at the mean bacterial spacer acquisition time vs $η / μ$ .

Points outlined in red are simulations where the ratio of large phage clones to bacterial clones exceeds 1.2. The phage fitness is the average per capita growth rate of phage clones conditioned on survival. Phage clones experience clonal interference at low $η$ and high $μ$ .

Figure 19

Download asset Open asset

Probability of the first spacer acquisition happening at time $t$ for four simulations with different values of $η$ and $C_{0} = 10^{4}, μ = 10^{- 5}$ , and $e = 0.95$ .

The mean of each distribution is shown as a vertical dashed line.

Figure 20

Download asset Open asset

Mean phage clone size at time of first spacer acquisition for simulation data, the predicted with Equation 31, and the prediction with the mode of the distribution given by Equation 29 for four simulations with different values of $η$ and $C_{0} = 10^{4}, μ = 10^{- 5}$ , and $e = 0.95$ .

Figure 21

Download asset Open asset

Measured mean phage clone size at the time of first spacer acquisition vs. the prediction given by $e^{s_{0} ⟨ t ⟩}$ of Equation 31.

Error bars are the standard deviation across three or more independent simulations.

Figure 22

Download asset Open asset

Mean number of large phage clones vs. mean number of bacterial clones in simulations.

For each simulation, we take a subset of 15 evenly spaced timepoints at steady state and calculate the size and number of phage clones present. We scale the observed clone sized distribution with Equation 28 and calculate the mean number of large phage clones by multiplying the total number of clones with the fraction of large phage clones given by $\sum_{n} P_{n}^{large}$ . We use the simulation mean total population sizes to calculate s₀ and $δ_{0}$ in Equation 27. We obtain the mean number of bacterial clones by averaging the number of clones present at 15 evenly spaced timepoints at steady state. Error bars are the standard deviation across three or more independent simulations.

Figure 23

Download asset Open asset

Effective $e = \frac{e}{n_{B}^{s} n_{V}} \sum_{i} n_{V}^{i} n_{B}^{i}$ vs $e / m$ across all simulations where $m \geq 1$ on average.

Error bars are the standard deviation across three or more independent simulations. The solid black line is $y = x$ .

Figure 24

Download asset Open asset

Four simulations with $C_{0} = 10000$ , $μ = 10^{- 5}$ , and $e = 0.95$ .

From top to bottom, $η$ increases by a factor of 10 in each row, from $η = 10^{- 5}$ in the top row to $η = 10^{- 2}$ in the bottom row. The first two columns show clone size distributions combined from 15 time points between 2000 and 10,000 bacterial generations. Bacteria are in the left column and phages in the middle column. The third column shows the pairwise clone sizes of matching clones at the last sampled time point (9467 generations). The expected large phage clone size is the total phage population divided by the mean number of bacterial clones. The inverse overlap is $\frac{e}{e_{eff}} = \frac{n_{B}^{s} n_{V}}{\sum_{i} n_{V}^{i} n_{B}^{i}}$ , which we assume is $\approx m$ as shown in Figure 23. The dashed line in the third column indicates the line that clone size pairs would fall on if they were perfectly correlated, and the red star indicates the mean large clone size for phages and the mean clone size for bacteria.

Figure 25

Download asset Open asset

Measured simulation phage clone mutation rate, establishment fraction, and mean time to extinction as a function of the theoretical prediction for each.

Highlighted in grey are parameter combinations for which no theoretically predicted $m$ could be determined; the predicted quantity is instead calculated with the simulation mean $m$ . Error bar are the stnadard deviation across three or more independent simulations.

Figure 26

Download asset Open asset

Total phage, total bacteria, and fraction of bacteria with spacers as a function of average immunity (effective $e$ ).

Blue points are simulation results. Error bars are standard deviation across three or more independent simulations. The solid black line is the solution given by Equations 42–44 (from Bonsma-Fisher et al., 2018) with the parameter $e$ replaced by effective $e$ . The horizontal grey dashed line corresponds to the no-CRISPR ( $e = 0$ ) mean-field solution (derived in Bonsma-Fisher et al., 2018).

Figure 27

Download asset Open asset

$n_{V}$ , $n_{B}$ , and $ν$ vs. $η$ for $e = 0.05$ , approximating $ν \approx - d / c$ .

Figure 28

Download asset Open asset

$n_{V}$ , $n_{B}$ , and $ν$ vs $η$ for $e = 0.15$ , approximating $ν \approx \frac{- c + \sqrt{c^{2} - 4 b d}}{2 b}$ .

Figure 29

Download asset Open asset

$n_{V}$ , $n_{B}$ , and $ν$ vs. $η$ for $e = 0.8$ , approximating $ν \approx \frac{- b + \sqrt{b^{2} - 4 a c}}{2 a}$ .

Figure 30

Download asset Open asset

Fraction of bacteria with spacers $ν$ vs. effective $e$ .

The solid line is the full numerical theoretical solution. The dashed black line is given by Equation 58. Error bars are the standard deviation across three or more independent simulations.

Figure 31

Download asset Open asset

Population standard deviation divided by population mean as a function of the distance from the average immunity critical point given by Equation 58.

Insets show a smaller x-axis range for the same quantities. Total phage (top left), total bacteria (top right), fraction of bacteria with spacers (bottom left) and total nutrients (bottom right) are plotted. X error bars are the standard deviation across three or more independent simulations.

Figure 32

Download asset Open asset

Approximate phage time to extinction vs. numerically calculated theoretical time to extinction, where the approximate time to extinction is given by (Koskella, 2014).

The full theoretical predicted value of $m$ is used in the approximate expression.

Figure 33

Download asset Open asset

Approximate phage mutation rate for $e = 0$ vs theoretical phage mutation rate.

Figure 34

Download asset Open asset

Measured mean $m$ in simulations vs. $a$ as given by Equation 81.

The solid line is $m$ vs. $a$ solved numerically using Equation 113. The right panel shows the same but with the two lowest values of $η$ removed. Error bars are the standard deviation across three or more independent simulations.

Figure 35

Download asset Open asset

Phage and bacteria centre of mass distance from the original phage and bacteria sequences.

The centre of mass distance is plotted in blue for phage (left) and bacteria (centre). Grey circles represent the size of clonal subpopulations at each distance from the ancestor sequence (arbitrary scale). The third panel shows the weighted average distance of the population from the centre of mass at that timepoint, a measure of the spread in sequences present at any time. In this simulation $C_{0} = 10^{4}$ , $η = 10^{- 5}$ , $μ = 10^{- 6}$ , $e = 0.95$ , and mean $m = 2.8$ .

Figure 36

Download asset Open asset

Figure 37

Download asset Open asset

Phage and bacteria centre of mass distance from the centre of mass at time $t - Δ t$ (left) and the distance divided by the time interval $Δ t$ (right).

Distances are averaged over the entire simulation at steady state; error bars are standard deviation. Simulation parameters are $C_{0} = 10^{4}$ , $η = 10^{- 2}$ , $μ = 10^{- 6}$ , $e = 0.95$ , and mean $m = 14.5$ .

Figure 38

Download asset Open asset

Maximum distance from ancestor population for bacteria vs. phage.

The maximum distance is highly correlated, indicating that the bacteria population tracks the phage population closely. Error bars are the standard deviation across three or more independent simulations.

Figure 39

Download asset Open asset

Phage mutational distance per generation vs. initial phage mutant fitness for simulations with $e = 0.95$ , $η = 10^{- 3}$ , and $μ = 3 \times 10^{- 7}$ .

Error bars are the standard deviation across multiple independent simulations and are shown in the positive direction only.

Figure 40

Download asset Open asset

Average population distance from the centre of mass at steady state for bacteria (left) and phages (right) vs. mean $m$ for simulations with one original phage clone ancestor.

The dashed line is $\sqrt{m} - 1$ , a purely phenomenological choice. Error bars are the standard deviation across three or more independent simulations.

Figure 41

Download asset Open asset

Average population distance from the centre of mass for bacteria (left) and phages (right) vs. mean $m$ for simulations with 50 original phage clones.

The dashed line is $\sqrt{m} - 1$ and the solid line is $m - 1$ . Error bars are the standard deviation across three or more independent simulations.

Figure 42

Download asset Open asset

A frame from a simulation movie at 5000 generations with $C_{0} = 10^{4}$ , $η = 10^{- 3}$ , $μ = 10^{- 5}$ , $e = 0.95$ , initial $m = 1$ .

Phages are on the left, bacteria with spacers on the right.

Figure 43

Download asset Open asset

Dendrogram resulting from agglomerative clustering with the L1 norm and linking clusters using the minimum distance between members.

The number of clusters is determined with a cutoff at a distance of 2. $C_{0} = 10^{4}$ , $η = 10^{- 3}$ , $μ = 10^{- 5}$ , $e = 0.95$ , initial $m = 1$ .

Figure 44

Download asset Open asset

Figure 45

Download asset Open asset

Clan number and size over time in a simulation with $C_{0} = 10^{4}$ , $η = 10^{- 3}$ , $μ = 10^{- 5}$ , $e = 0.95$ , initial $m = 1$ .

Figure 46

Download asset Open asset

Clan number and size over time in a simulation with $C_{0} = 10^{4}$ , $η = 10^{- 3}$ , $μ = 10^{- 5}$ , $e = 0.95$ , initial $m = 10$ .

Figure 47

Download asset Open asset

Average clan number vs. $m$ for all simulations that begin with 10 clones.

The dashed lines are $m$ divided by the mean bacterial clan size ( $\approx 1.3$ ) and $m / 2$ . Error bars are the standard deviation across three or more independent simulations.

Figure 48

Download asset Open asset

Average clan number vs. $m$ for all simulations that begin with one clone with $μ \geq 10^{- 6}$ .

The dashed lines are $m$ divided by the mean bacterial clan size and $m / 2$ . Error bars are the standard deviation across three or more independent simulations.

Figure 49

Download asset Open asset

Phage infection success probability $p_{V}$ as a function of mutational distance between spacers and protospacers for different definitions and degrees of cross-reactivity.

Definitions are plotted for $p_{V} = 0.02$ , $e = 0.95$ .

Figure 50

Download asset Open asset

Average immunity vs. diversity with different degrees of cross-reactivity for simulations with $e = 0.95$ , $μ = 10^{- 6}$ .

Dashed lines are simulations with cross-reactivity, solid line is simulations without cross-reactivity. Error bars are the standard deviation across three or more independent simulations.

Figure 51

Download asset Open asset

Mean phage clone size (top) and mean bacteria clone size (bottom) relative to time of phage mutation for different definitions of $p_{V}$ .

These simulations begin with one initial phage clone and parameters $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 0.0001$ , $e = 0.95$ .

Figure 52

Download asset Open asset

Number of bacterial clones (top) and average immunity (bottom) for simulations beginning with either 1 phage clone (left) or 10 phage clones (right).

These simulations have parameters $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 0.0001$ , $e = 0.95$ .

Figure 53

Download asset Open asset

Phage clone size (top) and bacteria clone size (bottom) in a simulation with cross-reactivity (step function CRISPR effectiveness with $θ = 1$ ).

This simulation has one initial phage clone and parameters $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 10^{- 4}$ , and $e = 0.95$ . Later times in this simulation are shown in Figure 55, and earlier times in this simulation are shown in Figure 54.

Figure 54

Download asset Open asset

Figure 55

Download asset Open asset

Phage clone size (top) and bacteria clone size (bottom) in a simulation with cross-reactivity (step function CRISPR effectiveness with $θ = 1$ ) showing the switch between a travelling wave and low turnover regime.

This simulation has one initial phage clone and parameters $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 10^{- 4}$ , and $e = 0.95$ . Earlier times in this simulation are shown in Figure 54.

Figure 56

Download asset Open asset

Phage clone size (top) and bacteria clone size (bottom) for a short time window of the simulation shown in Figure 54.

Large phage and bacteria clones are numbered in the legend; these numbers correspond to the numbers in Figure 57.

Figure 57

Download asset Open asset

Phage infection success probability matrix for each clone shown in Figure 56.

Dark blue is low infection success, light blue is high infection success.

Figure 58

Download asset Open asset

Mean phage clone size (top) and mean bacteria clone size (bottom) relative to the time of phage clone mutation, either normalized to surviving clones (left) or averaged over all trajectories (right) for the simulation shown in Figure 54 ( $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 10^{- 4}$ ).

Clones in the travelling wave regime (4000–6200 generations, orange) grow much more quickly than clones in the initial low-turnover regime (1000–3200 generations, blue) or the final low turnover regime (7000–10,000 generations, green). The black dashed line is the mean clone size for a simulation with the same parameters but without cross-reactivity.

Figure 59

Download asset Open asset

A slow-switching cross-reactivity regime for the simulation shown in Figure 54.

(Left) A subset of clone trajectories for phages (top) and bacteria (bottom) in a simulation with cross-reactivity (step function CRISPR effectiveness with $θ = 1$ ). Three trajectories are highlighted and coloured to show increasing time. The population is in a regime where matching clones are offset: because clones that are one mutation apart have the same complete immune overlap, the highlighted clones have large bacterial clone size well after the matching phage clone goes extinct. (Right) The clone size of the highlighted trajectories shown as a function of the bacterial marginal immunity against a particular phage clone (top) or the marginal immunity of a particular bacteria clone (bottom).

Figure 60

Download asset Open asset

Cross-reactivity leads to persistent oscillations.

(A) A subset of clone trajectories for phages (top) and bacteria (bottom) in a simulation with no cross-reactivity. Transient oscillations occur. One trajectory is highlighted and coloured to show increasing time. (B) The highlighted trajectory in (A) is shown as a function of the marginal immunity for phages (top) and bacteria (bottom). Clones experience an oscillating fitness that depends on their overlap from the other population. Arrows indicate the direction of increasing time in the oscillation. (C) A subset of clone trajectories for phages (top) and bacteria (bottom) in a simulation with cross-reactivity (step function CRISPR effectiveness with $θ = 1$ ). The population is in a regime where several clones experience persistent and rapid oscillations. One trajectory is highlighted and coloured to show increasing time. (D) The highlighted trajectory in (C) is shown as a function of the marginal immunity for phages (top) and bacteria (bottom). Clones experience an oscillating fitness that depends on their overlap from the other population. Arrows indicate the direction of increasing time in the oscillation. For all simulations $C_{0} = 10^{4}, e = 0.95, η = 10^{- 4}$ , and $μ = 10^{- 5}$ .

Figure 61

Download asset Open asset

Phage clone phylogenies for four simulations with different cross-reactivities and a higher mutation rate than shown in main text Figure 4: no cross-reactivity (A) and step-function cross-reactivity with $θ = 1$ (B, $θ = 2$ ) (C), and $θ = 3$ (D).

All simulations share all other parameters: $C_{0} = 10^{4}, η = 10^{- 4}, μ = 10^{- 5}, e = 0.95$ . Phage clones are plotted at the first time they pass a population size of 2 (to remove clutter from many new mutations destined for extinction), and the size of each circle is logarithmically proportional to the maximum size reached by that clone. Colours indicate the time of extinction of each clone. For each simulation with cross-reactivity, the left inset shows phage (top) and bacteria (bottom) clone sizes over time; colours indicate unique clone identities.

Figure 62

Download asset Open asset

Phage clone size (top) and bacteria clone size (bottom) for a short time window of the $θ = 2$ simulation shown in Figure 61C.

Each of the large trajectories oscillating out of phase between 5000 and 7000 generations is at least three mutations away from all of the others; they are all outside each other’s cross-reactivity radius.

Figure 63

Download asset Open asset

Matrix of mutational distance between each of the four largest phage clones shown in Figure 62; colours of those trajectories are labelled on the y axis.

Each clone is at least three mutations away from all other large clones.

Figure 64

Download asset Open asset

Time-shifted average immunity for three regimes of the simulation shown in Figure 54 ( $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 10^{- 4}$ , and $e = 0.95$ , step-function cross-reactivity with $θ = 1$ ).

The initial low-diversity regime (1000–3200 generations) and the low turnover, high diversity regime (7000–10,000 generations) had extremely low turnover, while the travelling wave regime (4000–6200 generations) had high average immunity near 0 delay that rapidly decayed to zero both in the past and future.

Figure 65

Download asset Open asset

Time-shifted average immunity for the corresponding simulation to Figure 64 without cross-reactivity ( $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 10^{- 4}$ , and $e = 0.95$ ). Peak average immunity is low because of high diversity, and immunity decays very gradually to zero in both the past and future.

Figure 66

Download asset Open asset

Time-shifted average immunity for four simulations with the same parameters but different types of cross-reactivity: no cross-reactivity (top), exponential cross-reactivity (middle top), step-function cross-reactivity with $θ = 1$ (middle bottom), and $θ = 2$ (bottom).

Shared parameters are $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 10^{- 4}$ , and $e = 0.95$ . Only the travelling-wave regime of each simulation with cross-reactivity was used to compare turnover in this regime. The grey shaded region is the standard deviation across averaged data.

Figure 67

Download asset Open asset

Average clan size across simulations with different parameters for different degrees of cross-reactivity with $m_{init} = 10$ .

Error bars are the standard deviation across three or more independent simulations.

Figure 68

Download asset Open asset

Average immunity vs. bacterial spacer array length for simulated distributions of protospacers and spacers.

We simulate 50 phages, each with 5 protospacers represented by letters from the alphabet, uniformly sampled. We simulate 420 bacterial spacers, drawn from the alphabet either uniformly (top row) or following an exponential distribution with mean 6 (bottom row). We construct arrays by sampling from the set of 420 spacers without replacement, either creating arrays of constant length (left column), or of variable length with a gaussian distribution (middle column) or exponential distribution (right column). Average immunity is calculated as in Equation 90 except the indices run over all arrays and not over individual sequences, and the presence of any matching pair gives perfect immunity. Blue points are average results over 50 simulations; error bars are standard deviation. The black dashed curve is given by $a_{n} = 1 - C^{n} = 1 - {(1 - a_{1})}^{n}$ .

Figure 69

Download asset Open asset

Average immunity decreases with coupled CRISPR diversity.

(A) Inferred initial average immunity as a function of diversity in the experiment of Common et al., 2020. Bacterial diversity was manipulated in the experiment by combining different numbers of bacterial clones defined by unique CRISPR spacers with a single phage strain that was able to infect only one of the clones. We calculated the expected initial average immunity based on the number of clones. (**B–D**) Results of a toy model of sorting a set of protospacers and spacers into arrays of variable length, either assigning protospacers fully randomly (B), as sequential mutations (C), or changing the last protospacer in the array only (D).

Figure 70

Download asset Open asset

Average immunity vs. diversity (number of unique types) for simulated distributions of spacers and protospacers with different bacterial array sizes (increasing top to bottom) and different numbers of protospacers per phage (increasing left to right).

Bacterial array sizes are drawn from a Gaussian distribution with mean given by the array mean for each row and a standard deviation of 2. Points are averages across 50 independent runs.

Figure 71

Download asset Open asset

Unique spacer types detected in our analysis of data from Paez-Espino et al., 2015 after grouping by 85% average similarity and removing single spacer counts (blue bars).

Counts reported in Paez-Espino et al., 2015 are red points.

Figure 72

Download asset Open asset

Total reads per time point that match phage (blue) or bacteria (orange).

Top: total reads. Middle: total reads divided by phage and bacteria genome sizes. Bottom: fraction of total reads matching bacteria or phages.

Figure 73

Download asset Open asset

Number of shared spacer types between bacteria and phage (top left), ratio of the number of phage types to the number of bacteria types (top right), number of bacteria types (bottom left), and number of phage types (bottom right) as a function of the sampling date for data we analysed from Paez-Espino et al., 2015.

Colours indicate different similarity grouping thresholds. Spacer counts include wild-type spacers, and protospacers are included only if they possess a perfect PAM. Data is summed over CRISPR1 and CRISPR3. All sequences are included regardless of abundance.

Figure 74

Download asset Open asset

Total phage and bacteria population size for the MOI2B experiment in Paez-Espino et al.Paez-Espino et al., 2015.

Circles are points digitized from Figure 1A in Paez-Espino et al., 2015; squares are the population size interpolated to match the sequencing dates in the experiment.

Figure 75

Download asset Open asset

Total number of unique protospacer sequences after removing all sequences that cluster with wild-type sequences at different similarity thresholds.

Only protospacers with a perfect PAM are included.

Figure 76

Download asset Open asset

Phage (top) and bacteria (bottom) population size vs. average immunity calculated with data from Paez-Espino et al., 2015.

Protospacers are included if they have a perfect PAM, and all wild-type spacers are included. Spacers and protospacers are grouped with an 85% similarity threshold. Colours from blue to yellow indicate increasing time.

Figure 77

Download asset Open asset

Average immunity vs time in the MOI-2B experiment.

Protospacers are included if they have a perfect PAM, and all wild-type spacers are included. Spacers and protospacers are grouped with an 85% similarity threshold.

Figure 78

Download asset Open asset

Histogram of reads that matched both the *Gordonia* MAG and either phage reference genome vs the base-10 log difference in e-value between the matches for accession SRR9260993.

A positive value means the bacteria match has a lower e-value than the phage match. The vertical dashed line indicates the -25 cutoff; matches to the left were considered phage matches, matches to the right were considered bacteria matches.

Figure 79

Download asset Open asset

Probability logo for the four nucleotides at the 5’ end of potential protospacer sequences, generated with WebLogo; the spacer sequence starts at the right edge of the logo meaning that the reverse-complement PAM is GTT.

Figure 80

Download asset Open asset

Phage genome (top) and bacteria genome (bottom) coverage, digitized from Figure 2D of Guerrero et al., 2021a vs. number of reads assigned to phage or bacteria from our analysis.

Each marker is a separate time point.

Figure 81

Download asset Open asset

Abundance over time for the 20 largest bacteria clones (top) and phage clones (bottom) over time for data from Guerrero et al., 2021a.

The left panels show absolute counts, and the right panels show fractional abundance for the included types. None of the top types are shared between phage and bacteria.

Figure 82

Download asset Open asset

Figure 83

Download asset Open asset

Fraction of spacer types (left) and protospacer types (right) remaining as a function of time delay, averaged over the entire time series from Guerrero et al., 2021a.

Types are grouped with an 85% similarity threshold. The grey shaded region is the standard deviation across averaged data.

Figure 84

Download asset Open asset

Average immunity at each time point in the time series.

Raw values cannot be larger than 1, but plotted values are multiplied by 1956, the average number of protospacers with the GTT PAM from the phage DC-56 and DS-92 genomes, yielding some values larger than 1.

Figure 85

Download asset Open asset

Distribution of average immunity values at zero time shift (blue) and at 500 days time shift for future phages (orange) and 500 days time shift for past phages (green).

Figure 86

Download asset Open asset

p-values for the Wilcoxon signed-rank test comparing the average immunity at zero time shift with all other time shifts.

Bacteria immunity against past phages is generally not significantly lower than bacterial immunity against current phages (blue), while bacterial immunity against future phages is significantly lower than immunity against current phages for almost all time shifts (orange). The time shifts are not symmetric; bacterial overlap with past phages is not necessarily the same as phage overlap with future bacteria and vice versa. The number of points that are available to compare decreases as the time shifts get larger (black dashed line).

Figure 87

Download asset Open asset

Time-shifted average immunity for each time point as a function of the zero-shift average immunity of that data point for shifts comparing past phages (top row) and future phages (bottom row).

Figure 88

Download asset Open asset

Phage (top) and bacteria (bottom) population size vs. average immunity calculated with data from Guerrero et al., 2021a.

Protospacers are included if they have a perfect PAM, and all wild-type spacers are included. Spacers and protospacers are grouped with an 85% similarity threshold. Colours from blue to yellow indicate increasing time.

Appendix 2—figure 1

Download asset Open asset

Diversity vs. $η$ (left), $μ$ (centre), and $e$ (right) for $C_{0} = 10^{4}$ .

The $η$ dependence of diversity is not very well predicted even by the full numerical solution, but for mutation rate and spacer effectiveness the approximate solutions do pretty well in this regime. The simulation data increase in diversity as a function of spacer acquisition probability actually goes more like $m \propto \ln η$ .

Appendix 2—figure 2

Download asset Open asset

Approximate predictions for diversity in our coevolution model (blue) or a simple model with mutation but no coevolution (orange).

(A) Predicted diversity as a function of spacer acquisition probability in our coevolution model as given by Equation 111 for $C_{0} = 10^{4}$ , $e = 0.95$ , and $μ = 10^{- 6}$ (blue). Predicted diversity in the non-coevolving model as given by $m = - N_{b} \frac{u}{1 - u} \ln u$ with $N_{b} = 1000$ (orange). The low- $μ$ limit is $m = - N_{b} u \ln u$ . The high- $μ$ limit is a series expansion in $ϵ$ for $μ = 1 - ϵ$ giving $m \approx N_{b} (\frac{1}{3} + \frac{5 μ}{6} - \frac{μ^{2}}{6})$ . (B) Fold change in diversity as a function of fold change in mutation rate or spacer acquisition probability under the simple approximation that $m \sim η^{1 / 3}$ (blue) and $m \sim u$ (orange).

Appendix 2—figure 3

Download asset Open asset

Fold change in diversity (number of species) as a function of mutation rate and $k$ , the fold change in mutation rate.

(A) Fold change in diversity as a function of fold change in mutation rate in the coevolution model: diversity increases by approximately a factor of $k^{\frac{1}{3}}$ , independent of mutation rate (solid line). (B) Fold change in diversity as a function of both mutation rate and fold change in mutation rate in the model without coevolution.

Appendix 3—figure 1

Download asset Open asset

Phage clone initial growth rate vs. total bacteria normalized by the initial nutrient concentration C₀.

Phage clone growth rate is as defined in Figure 3—figure supplement 1; for each simulation, the average phage clone growth rate is the derivative of the average phage clone size, averaged across all trajectories after steady state; plotted points and error bars are the average across three or more simulations. The phage initial fitness depends slightly on the phage mutation rate (mutants decrease the growth rate of a particular phage clone), but this dependence is slight enough that all mutation rates collapse onto the theoretical line. Here we plot data with $μ = 10^{- 6}$ . The phage clone initial growth rate also does not depend on $e$ or $η$ because new phage mutants see the bacteria population as if it did not have any CRISPR immunity. The theoretical initial phage clone growth rate is given by Equation 122. The effective lower bound of $n_{B} / C_{0}$ is set by the steady-state population size without CRISPR immunity: $n_{B} / C_{0} = \frac{f g}{α (B p_{V} - 1)} \approx 0.15$ .

Appendix 3—figure 2

Download asset Open asset

Clone size histograms (left) and cumulative distributions (right) for four different values of the spacer acquisition probability $η$ .

In all simulations, $C_{0} = 10^{4}$ , $e = 0.95$ , and $μ = 3 \times 10^{- 6}$ . We sample 30 evenly spaced time points between 2000 and 10,000 bacterial generations and combine the clone sizes at each of the sampled points to create the clone size distributions plotted. Solid lines show three different theoretical solutions. The solid blue line is given by Equation 133 with all population quantities predicted from solving the system of Equations 13–17 with $m$ given by Equation 32 and $n_{V}^{i} = n_{V} / m$ . The solid orange line is given by Equation 133, with the value $n_{V}^{i}$ determined by maximum likelihood estimation (MLE) to give the best fit to the data. For the two largest values of $η$ , the value of $n_{V}^{i}$ returned by the MLE fit is smaller than the theoretical value of $n_{V}^{i}^{*}$ , while for the two smallest values of $η$ the MLE value of $n_{V}^{i}$ is larger. For large enough values of $n_{V}^{i}$ , the bacteria clone death rate $d$ is smaller than the birth rate $b$ which violates the assumptions used to derive Equation 133. This happens for the MLE fit at the two smallest $η$ values and hence no MLE solution is plotted. The solid green line is an average distribution calculated by solving Equation 133 at each observed large phage clone size and averaging across the distribution of clone sizes; i.e. $P (n_{B}^{i}) = P (n_{B}^{i} | n_{V}^{i}) P (n_{V}^{i})$ . The large phage clone distribution is given by Equation 28.

Appendix 3—figure 3

Download asset Open asset

Bacteria and phage clone trajectories aligned to the time at which bacteria trajectories go extinct.

Bacteria trajectories are included if they reach size $n_{B}^{i}^{*}$ given by Equation 25 and all corresponding phage trajectories are plotted. In this simulation $C_{0} = 3 \times 10^{4}$ , $e = 0.8$ , $η = 10^{- 3}$ , and $μ = 10^{- 6}$ .

Appendix 3—figure 4

Download asset Open asset

Mean time to extinction for bacterial clones after reaching size $n_{B}^{i}^{*}$ as a function of $η$ for $C_{0} = 10^{4}$ , $e = 0.95$ , and $μ = 3 \times 10^{- 6}$ .

The solid line is given by Equation 139 with $n = n_{B}^{i}^{*}$ , and the dashed line is given by numerically solving Equation 137.

Appendix 3—figure 5

Download asset Open asset

Predicted bacteria clone time to extinction with drift.

Measured vs. predicted mean time to extinction for bacterial clones after reaching size $n_{B^{i}}^{*}$ . The predicted time to extinction is the solution with drift, given by numerically solving Equation 167.

Appendix 3—figure 6

Download asset Open asset

Approximations for bacteria clone extinction.

Measured mean time to extinction for large bacterial clones vs. three successively more aggressive analytic approximations for the mean time to extinction. (A) The predicted time to extinction is given by taking a series expansion in $e / m$ and keeping terms to 1st order. (B) A series expansion in $e / m$ and keeping the 0th order term ( $e = 0$ ) only. (C) Series expansion in $e / m$ and $r$ and keeping the 0th order term only (Equation 143).

Appendix 3—figure 7

Download asset Open asset

Approximations for bacteria clone extinction as a function of mean clone size.

Bacteria extinction time as a function of mean clone size. Large bacteria clone mean time to extinction as a function of $\frac{n_{B}}{m}$ , where $m$ is measured from simulations and $n_{B}$ is the total bacteria population size. This parameter combination describes the trend in both phage and bacteria extinction reasonably well at large values of (bottom panels). The solid line is given by Equation 145 without the $1 + \ln m$ term.

Appendix 3—figure 8

Download asset Open asset

Phage clone extinction times and theoretical predictions for a simulation with parameters $C_{0} = 10^{4}, η = 0.001, e = 0.95, μ = 10^{- 5}$ .

Time zero for each trajectory is the time at which that clone arose by mutation. All simulation trajectories are plotted in blue, and a subset of trajectories that do not reach a size of $n_{V}^{i}^{*} = 16170$ as given by Equation 26 are plotted in orange. Trajectories that do not become established by this definition go extinct more quickly than ones that do. All other curves are theoretical predictions for the extinction time. The green and red solid lines and purple and brown dashed lines show a numerical solution to Equation 153 with different values for $s$ . The remaining dashed lines are a small time approximation given by Equation 163, and a large time approximation given by Equation 156 with either $s = s_{0}$ or $s = 0$ . All of these predictions agree with the simulation data in different regimes; none accurately captures the entire timecourse of extinctions.

Appendix 3—figure 9

Download asset Open asset

Large clone extinction times from a simulation with parameters $C_{0} = 10^{4}, η = 0.001, e = 0.95, μ = 10^{- 5}$ .

Trajectories are counted as large if the phage clone size passes $n_{V}^{i}^{*} = 16170$ , the theoretical deterministic mean phage clone size for these parameters given by Equation 26. Once a trajectory reaches $n_{V}^{i}^{*}$ , we count that point as time zero to measure the extinction time of large clones. The numerical solution is given by solving Equation 153 with $s = 0$ ( $p = 1 / B$ ). The large time approximation is given by Equation 157. The orange dashed line gives an exponential decay prediction with the mean time to extinction given by Equation 171 with $n = n_{V}^{i}^{*}$ .

Appendix 3—figure 10

Download asset Open asset

Mean time to extinction for large bacteria (left) and phage (right) clones as a function of the neutral fitness mean time to extinction prediction.

The solid black lines describe approximate analytic expressions for bacteria and phage time to extinction using a neutral fitness assumption. For bacteria, the black line solves Equation 139 using the mean bacteria and phage clone sizes. For phages, the black line solves Equation 171 using the mean phage clone size in place of $n$ .

Appendix 3—figure 11

Download asset Open asset

Phage extinction time as a function of mean clone size.

Large phage clone mean time to extinction as a function of $\frac{n_{V}}{m}$ , where $m$ is measured from simulations and $n_{V}$ is the total phage population size. The solid line is given by Equation 172 without the $1 + \ln m$ term.

Tables

Table 1

Model reactions.

$b^{0, i} + C \overset{g}{\to} 2 b^{0, i}$	Bacterium divides
$b^{0, i} \overset{F}{\to} \emptyset$	Bacterium flows out
$V^{j} \overset{F}{\to} \emptyset$	Phage flows out
$\emptyset \overset{F C_{0}}{\to} C$	Nutrients flow in
$C \overset{F}{\to} \emptyset$	Nutrients flow out
$\sum_{n = 0}^{B} (b^{0} + V^{j} \overset{α p_{V} P_{n}}{\to} (B - n) V^{j} + \sum_{k = 1}^{n} V^{m + k})$	Interaction, phage wins,
	$P_{n}$ probability of $n$ mutuant phages
$b^{0} + V^{j} \to α (1 - p_{V}) (1 - η) b^{0}$	Interaction, bacterium survives
$b^{0} + V^{j} \to α (1 - p_{V}) η b^{i}$	Interaction, bacterium survives and acquires a spacer
$\sum_{n = 0}^{B} (b^{i} + V^{j} \to α p_{V} (i, j) P_{n} (B - n) V^{j} + \sum_{k = 1}^{n} V^{m + k})$	Interaction, phage wins,
	$P_{n}$ probability of $n$ mutant phages
$b^{i} + V^{j} \underset{\to}{α (1 - p_{V} (i, j))} b^{i}$	Interaction, bacterium survives
$b^{i} \overset{r}{\to} b^{0}$	Bacterium loses spacer

Table 2

Simulation length.

C₀	Simulation length (bacterial generations)	Steady-state start time ( $t_{s s}$ )
300	10,000	2000
1000	10,000	2000
3000	10,000	2000
10,000	10,000	2000
30,000	15,000	3000
100,000	20,000	4000
300,000	25,000	5000
1,000,000	30,000	6000

Table 3

Model parameters.

Parameter	Description	Value
$\frac{1}{g C_{0}}$	Bacterial doubling time	41.7 min
C₀	Inflow nutrient concentration in	10² to 10^-6
	Units of bacterial cell density
$α$	Phage adsorption rate	$2 \times 10^{- 2} / C_{0}$
$B$	Phage burst size	170
$F$	Chemostat flow rate	$0.3 g C_{0}$
$p_{V}$	Probability of phage success	0.02
	for bacteria without spacers
$e$	Spacer effectiveness	0.1 to 0.95
$r$	Rate of spacer loss	$0.04 g C_{0}$
$η$	Probability of spacer acquisition	10^-5 to 10^-2
$μ$	Phage mutation rate per base per generation	10^-8 to 10^-4
$L$	Phage protospacer length in nucleotides	30

Parameter values are as above unless otherwise indicated. Representative values estimated for Streptococcus thermophilus bacteria in lab conditions.

Table 4

$ν$ approximations.

$η$	$e_{effective}$
	≤0.01	0.01 to 0.05	0.05 to 0.1	0.1 to 0.5	0.5 to 1

$\leq 10^{- 5}$	$- \frac{d}{c}$	$- \frac{d}{c}$	$- \frac{d}{c}$	$\frac{- b + \sqrt{b^{2} - 4 a c}}{2 a}$	$\frac{- b + \sqrt{b^{2} - 4 a c}}{2 a}$
10^-5 to10^-4	$- \frac{d}{c}$	$- \frac{d}{c}$	$\frac{- c + \sqrt{c^{2} - 4 b d}}{2 b}$	$\frac{- b + \sqrt{b^{2} - 4 a c}}{2 a}$
10^-4 to 10^-2	$- \frac{d}{c}$	$- \frac{d}{c}$	$\frac{- c + \sqrt{c^{2} - 4 b d}}{2 b}$
$\geq 10^{- 2}$	$- \frac{d}{c}$	$\frac{- c + \sqrt{c^{2} - 4 b d}}{2 b}$	$\frac{- c + \sqrt{c^{2} - 4 b d}}{2 b}$

Note: the drop- $a$ solution also works for below $e = 0.1$ for all values of $η$ (since ${\frac{2 b d}{c}}^{2} \approx 0$ ), but the $- d / c$ solution is simpler and so preferred for the very low $e$ range.

Additional files

MDAR checklist: https://cdn.elifesciences.org/articles/81692/elife-81692-mdarchecklist1-v2.docx
Download elife-81692-mdarchecklist1-v2.docx

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Madeleine Bonsma-Fisher
Sidhartha Goyal

(2023)

Dynamics of immune memory and learning in bacterial communities

eLife 12:e81692.

https://doi.org/10.7554/eLife.81692

Share this article

Cite this article

Model description.

Probability of stochastic extinction at low spacer acquisition.

Probability of stochastic extinction at high spacer acquisition.

Time to extinction for phages vs. mutation rate.

Time to extinction for bacteria vs. mutation rate.

Time to extinction for phages for different initial diversity and low spacer acquisition.

Time to extinction for phages for different initial diversity and high spacer acquisition.

An animation of a typical simulation of bacteria and phages interacting with CRISPR immunity.

Diversity depends sub-linearly on parameters.

Predicted diversity grouped by spacer effectiveness.

Predicted diversity grouped by mutation rate.

Clone size histograms by spacer effectiveness.

Clone size histograms by phage mutation rate.

The fate of individual clones.

Average clone sizes and fitness over time in a simulation.

Approximations for phage establishment probability.

Theoretical phage establishment probability with approximations.

Cross-reactivity leads to ‘spindly’ phylogenies and regime switching.

Animation of simulation without cross-reactivity for Figure 4A.

Simulation animation with exponential cross-reactivity for Figure 4B.

Simulation animation with step-function cross-reactivity for Figure 4C.

Simulation animation with step-function cross-reactivity for Figure 4D.

Average immunity underlies population outcomes.

Total population sizes predicted by average immunity in a simulation with cross-reactivity.

Phage evolution and spacer turnover.

Clone size PCA for a simulation with exponential cross-reactivity.

Genetic distance vs. diversity and time to extinction.

Spacer turnover by spacer effectiveness.

Spacer turnover by phage mutation rate.

PCA decomposition of phage and bacteria protospacer and spacer clone abundances for a simulation of bacteria and phages interacting with CRISPR immunity.

Quantifying immune memory in data.

Time-shifted average immunity for experimental coevolution data.

Time-shifted average immunity for experimental coevolution data with partial PAM matches included.

Time-shifted average immunity for experimental coevolution data with all sequences regardless of PAM.

Time-shifted average immunity for experimental coevolution data with data trimmed from start and end.

Time-shifted average immunity for wastewater with data trimmed from start and end.

Bootstrapped average immunity with all points randomly shuffled.

Bootstrapped average immunity with pairs of bacteria and phage clone sizes shuffled.

Time-shifted average immunity for wastewater with data trimmed from the end.

Time-shifted average immunity for wastewater with data trimmed from the start.

Time-shifted average immunity for wastewater with data trimmed from the start and end.

Time-shifted average immunity by spacer effectiveness.

Time-shifted average immunity by phage mutation rate.

Total phage, total bacteria (left), and mean number of bacterial clones (right) vs. simulation time for five simulations with μ=10-6, e=0.5, η=0.001, and C0 ranging from 300 (top row) to 30,000 (bottom row).

Mean number of bacterial clones after t=ts⁢s bacterial generations vs. initial number of phage clones for e=0.8, η=10-3.

Phage (A) and bacteria (B) mean clone size in a simulation, either conditioned on survival (blue circles) or including extinct clones (orange circles).

Total phage (nV) and total bacteria (nB) as a function of time for a Gillespie simulation and a tau-leaping simulation with the same parameters.

Number of bacterial clones (m) vs. simulation time for three sets of simulations, each beginning with 1, 10, or 50 phage clones.

Average population size of phage clones (solid lines) and bacterial clones (dashed lines) in a Gillespie and tau-leaping simulation with the same parameters.

10 spacer trajectories for a Gillespie simulation (A) and tau-leaping simulation (B).

Phage clone size distribution from 15 combined time points for a simulation with the parameters C0=104, e=0.95, η=0.01, and μ=3×10-6.

Mean phage clone size at the time of first spacer acquisition vs. the deterministic mean phage clone size.

Mean time to extinction for phage clones vs. the timescale of bacteria spacer acquisition given by 1/D where D=α⁢η⁢(1-pV)⁢nVi*⁢nB0⁢g⁢C0.

Ratio of average phage clone initial fitness to the average phage clone fitness at the mean bacterial spacer acquisition time vs η/μ.

Probability of the first spacer acquisition happening at time t for four simulations with different values of η and C0=104,μ=10-5, and e=0.95.

Mean phage clone size at time of first spacer acquisition for simulation data, the predicted with Equation 31, and the prediction with the mode of the distribution given by Equation 29 for four simulations with different values of η and C0=104,μ=10-5, and e=0.95.

Measured mean phage clone size at the time of first spacer acquisition vs. the prediction given by es0⟨t⟩ of Equation 31.

Mean number of large phage clones vs. mean number of bacterial clones in simulations.

Effective e=enBs⁢nV⁢∑inVi⁢nBi vs e/m across all simulations where m≥1 on average.

Four simulations with C0=10000, μ=10-5, and e=0.95.

Measured simulation phage clone mutation rate, establishment fraction, and mean time to extinction as a function of the theoretical prediction for each.

Total phage, total bacteria, and fraction of bacteria with spacers as a function of average immunity (effective e).

nV, nB, and ν vs. η for e=0.05, approximating ν≈-d/c.

nV, nB, and ν vs η for e=0.15, approximating ν≈-c+c2-4⁢b⁢d2⁢b.

nV, nB, and ν vs. η for e=0.8, approximating ν≈-b+b2-4⁢a⁢c2⁢a.

Fraction of bacteria with spacers ν vs. effective e.

Population standard deviation divided by population mean as a function of the distance from the average immunity critical point given by Equation 58.

Approximate phage time to extinction vs. numerically calculated theoretical time to extinction, where the approximate time to extinction is given by (Koskella, 2014).

Approximate phage mutation rate for e=0 vs theoretical phage mutation rate.

Measured mean m in simulations vs. a as given by Equation 81.

Phage and bacteria centre of mass distance from the original phage and bacteria sequences.

Phage and bacteria centre of mass distance from the original phage and bacteria sequences.

Phage and bacteria centre of mass distance from the centre of mass at time t-Δ⁢t (left) and the distance divided by the time interval Δ⁢t (right).

Maximum distance from ancestor population for bacteria vs. phage.

Phage mutational distance per generation vs. initial phage mutant fitness for simulations with e=0.95, η=10-3, and μ=3×10-7.

Average population distance from the centre of mass at steady state for bacteria (left) and phages (right) vs. mean m for simulations with one original phage clone ancestor.

Average population distance from the centre of mass for bacteria (left) and phages (right) vs. mean m for simulations with 50 original phage clones.

A frame from a simulation movie at 5000 generations with C0=104, η=10-3, μ=10-5, e=0.95, initial m=1.

Total phage, total bacteria (left), and mean number of bacterial clones (right) vs. simulation time for five simulations with $μ = 10^{- 6}$ , $e = 0.5$ , $η = 0.001$ , and C₀ ranging from 300 (top row) to 30,000 (bottom row).

Mean number of bacterial clones after $t = t_{s s}$ bacterial generations vs. initial number of phage clones for $e = 0.8$ , $η = 10^{- 3}$ .

Total phage ( $n_{V}$ ) and total bacteria ( $n_{B}$ ) as a function of time for a Gillespie simulation and a tau-leaping simulation with the same parameters.

Number of bacterial clones ( $m$ ) vs. simulation time for three sets of simulations, each beginning with 1, 10, or 50 phage clones.

Phage clone size distribution from 15 combined time points for a simulation with the parameters $C_{0} = 10^{4}$ , $e = 0.95$ , $η = 0.01$ , and $μ = 3 \times 10^{- 6}$ .

Mean time to extinction for phage clones vs. the timescale of bacteria spacer acquisition given by $1 / D$ where $D = α η (1 - p_{V}) n_{V}^{i}^{*} n_{B}^{0} g C_{0}$ .

Ratio of average phage clone initial fitness to the average phage clone fitness at the mean bacterial spacer acquisition time vs $η / μ$ .

Probability of the first spacer acquisition happening at time $t$ for four simulations with different values of $η$ and $C_{0} = 10^{4}, μ = 10^{- 5}$ , and $e = 0.95$ .

Mean phage clone size at time of first spacer acquisition for simulation data, the predicted with Equation 31, and the prediction with the mode of the distribution given by Equation 29 for four simulations with different values of $η$ and $C_{0} = 10^{4}, μ = 10^{- 5}$ , and $e = 0.95$ .

Measured mean phage clone size at the time of first spacer acquisition vs. the prediction given by $e^{s_{0} ⟨ t ⟩}$ of Equation 31.

Effective $e = \frac{e}{n_{B}^{s} n_{V}} \sum_{i} n_{V}^{i} n_{B}^{i}$ vs $e / m$ across all simulations where $m \geq 1$ on average.

Four simulations with $C_{0} = 10000$ , $μ = 10^{- 5}$ , and $e = 0.95$ .

Total phage, total bacteria, and fraction of bacteria with spacers as a function of average immunity (effective $e$ ).

$n_{V}$ , $n_{B}$ , and $ν$ vs. $η$ for $e = 0.05$ , approximating $ν \approx - d / c$ .

$n_{V}$ , $n_{B}$ , and $ν$ vs $η$ for $e = 0.15$ , approximating $ν \approx \frac{- c + \sqrt{c^{2} - 4 b d}}{2 b}$ .

$n_{V}$ , $n_{B}$ , and $ν$ vs. $η$ for $e = 0.8$ , approximating $ν \approx \frac{- b + \sqrt{b^{2} - 4 a c}}{2 a}$ .

Fraction of bacteria with spacers $ν$ vs. effective $e$ .

Approximate phage mutation rate for $e = 0$ vs theoretical phage mutation rate.

Measured mean $m$ in simulations vs. $a$ as given by Equation 81.

Phage and bacteria centre of mass distance from the centre of mass at time $t - Δ t$ (left) and the distance divided by the time interval $Δ t$ (right).

Phage mutational distance per generation vs. initial phage mutant fitness for simulations with $e = 0.95$ , $η = 10^{- 3}$ , and $μ = 3 \times 10^{- 7}$ .

Average population distance from the centre of mass at steady state for bacteria (left) and phages (right) vs. mean $m$ for simulations with one original phage clone ancestor.

Average population distance from the centre of mass for bacteria (left) and phages (right) vs. mean $m$ for simulations with 50 original phage clones.

A frame from a simulation movie at 5000 generations with $C_{0} = 10^{4}$ , $η = 10^{- 3}$ , $μ = 10^{- 5}$ , $e = 0.95$ , initial $m = 1$ .

Clan number and size over time in a simulation with $C_{0} = 10^{4}$ , $η = 10^{- 3}$ , $μ = 10^{- 5}$ , $e = 0.95$ , initial $m = 1$ .

Clan number and size over time in a simulation with $C_{0} = 10^{4}$ , $η = 10^{- 3}$ , $μ = 10^{- 5}$ , $e = 0.95$ , initial $m = 10$ .

Average clan number vs. $m$ for all simulations that begin with 10 clones.

Average clan number vs. $m$ for all simulations that begin with one clone with $μ \geq 10^{- 6}$ .

Phage infection success probability $p_{V}$ as a function of mutational distance between spacers and protospacers for different definitions and degrees of cross-reactivity.

Average immunity vs. diversity with different degrees of cross-reactivity for simulations with $e = 0.95$ , $μ = 10^{- 6}$ .

Mean phage clone size (top) and mean bacteria clone size (bottom) relative to time of phage mutation for different definitions of $p_{V}$ .

Phage clone size (top) and bacteria clone size (bottom) in a simulation with cross-reactivity (step function CRISPR effectiveness with $θ = 1$ ).

Phage clone size (top) and bacteria clone size (bottom) in a simulation with cross-reactivity (step function CRISPR effectiveness with $θ = 1$ ).

Phage clone size (top) and bacteria clone size (bottom) in a simulation with cross-reactivity (step function CRISPR effectiveness with $θ = 1$ ) showing the switch between a travelling wave and low turnover regime.

Mean phage clone size (top) and mean bacteria clone size (bottom) relative to the time of phage clone mutation, either normalized to surviving clones (left) or averaged over all trajectories (right) for the simulation shown in Figure 54 ( $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 10^{- 4}$ ).

Phage clone phylogenies for four simulations with different cross-reactivities and a higher mutation rate than shown in main text Figure 4: no cross-reactivity (A) and step-function cross-reactivity with $θ = 1$ (B, $θ = 2$ ) (C), and $θ = 3$ (D).

Phage clone size (top) and bacteria clone size (bottom) for a short time window of the $θ = 2$ simulation shown in Figure 61C.

Time-shifted average immunity for three regimes of the simulation shown in Figure 54 ( $C_{0} = 10^{4}$ , $μ = 10^{- 6}$ , $η = 10^{- 4}$ , and $e = 0.95$ , step-function cross-reactivity with $θ = 1$ ).

Time-shifted average immunity for four simulations with the same parameters but different types of cross-reactivity: no cross-reactivity (top), exponential cross-reactivity (middle top), step-function cross-reactivity with $θ = 1$ (middle bottom), and $θ = 2$ (bottom).

Average clan size across simulations with different parameters for different degrees of cross-reactivity with $m_{init} = 10$ .

Diversity vs. $η$ (left), $μ$ (centre), and $e$ (right) for $C_{0} = 10^{4}$ .

Fold change in diversity (number of species) as a function of mutation rate and $k$ , the fold change in mutation rate.

Phage clone initial growth rate vs. total bacteria normalized by the initial nutrient concentration C₀.

Clone size histograms (left) and cumulative distributions (right) for four different values of the spacer acquisition probability $η$ .

Mean time to extinction for bacterial clones after reaching size $n_{B}^{i}^{*}$ as a function of $η$ for $C_{0} = 10^{4}$ , $e = 0.95$ , and $μ = 3 \times 10^{- 6}$ .

Phage clone extinction times and theoretical predictions for a simulation with parameters $C_{0} = 10^{4}, η = 0.001, e = 0.95, μ = 10^{- 5}$ .

Large clone extinction times from a simulation with parameters $C_{0} = 10^{4}, η = 0.001, e = 0.95, μ = 10^{- 5}$ .

$ν$ approximations.