Dynamics of immune memory and learning in bacterial communities

  1. Madeleine Bonsma-Fisher
  2. Sidhartha Goyal  Is a corresponding author
  1. Department of Physics, University of Toronto, Canada
  2. Institute of Biomaterials & Biomedical Engineering, University of Toronto, Canada
102 figures, 4 tables and 1 additional file

Figures

Figure 1 with 7 supplements
Model description.

(A) We model bacteria and phages interacting in a well-mixed vessel. We track nutrient concentration, phage population size (nV), and bacteria population size (nB). Bacteria can either have no …

Figure 1—figure supplement 1
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 …

Figure 1—figure supplement 2
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 …

Figure 1—figure supplement 3
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 …

Figure 1—figure supplement 4
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 …

Figure 1—figure supplement 5
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 …

Figure 1—figure supplement 6
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 …

Figure 1—video 1
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 …

Figure 2 with 4 supplements
Diversity depends sub-linearly on parameters.

(A) Bacteria and phage clone size distributions normalized to the measured mean clone size for C0=105, μ=3×10-7, and e=0.95. As η increases, both clone size distributions become more sharply peaked. (B) The mean …

Figure 2—figure supplement 1
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 …

Figure 2—figure supplement 2
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 …

Figure 2—figure supplement 3
Clone size histograms by spacer effectiveness.

Bacteria and phage clone size distributions normalized to the measured mean clone size for C0=105, μ=3×10-7, and η=0.001.

Figure 2—figure supplement 4
Clone size histograms by phage mutation rate.

Bacteria and phage clone size distributions normalized to the measured mean clone size for C0=105, e=0.95, and η=0.001.

Figure 3 with 3 supplements
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 …

Figure 3—figure supplement 1
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 …

Figure 3—figure supplement 2
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 …

Figure 3—figure supplement 3
Theoretical phage establishment probability with approximations.

Theoretical phage establishment probability vs. e/m. Solid lines are the numerical theoretical solution. Light dashed lines (‘small eeff’) are the approximate solution given by Equation 54 (with the …

Figure 4 with 4 supplements
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: C0=104,η=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
Animation of simulation without cross-reactivity for Figure 4A.
Figure 4—video 2
Simulation animation with exponential cross-reactivity for Figure 4B.
Figure 4—video 3
Simulation animation with step-function cross-reactivity for Figure 4C.
Figure 4—video 4
Simulation animation with step-function cross-reactivity for Figure 4D.
Figure 5 with 1 supplement
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 …

Figure 5—figure supplement 1
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 …

Figure 6 with 5 supplements
Phage evolution and spacer turnover.

(A) Principal Component Analysis (PCA) decomposition of phage and bacteria clone abundances for a simulation with C0=104, e=0.95, η=10-4, and μ=10-5. Clone abundances are normalized at each time point, then PCA is …

Figure 6—figure supplement 1
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 C0=104, e=0.95,,η=10-4 and.μ=10-6 Clone abundances are normalized at each time point, then PCA is …

Figure 6—figure supplement 2
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 …

Figure 6—figure supplement 3
Spacer turnover by spacer effectiveness.

Spacer turnover as a function of time delay for four simulations with C0=104, eta=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 …

Figure 6—figure supplement 4
Spacer turnover by phage mutation rate.

Spacer turnover as a function of time delay for four simulations with C0=104, e=0.95, and eta=0.001. The fraction of bacterial clones remaining is the fraction of clones that were present at time t that are still …

Figure 6—video 1
PCA decomposition of phage and bacteria protospacer and spacer clone abundances for a simulation of bacteria and phages interacting with CRISPR immunity.

Clone abundances are normalized at each time point, then PCA is performed for the entire phage time series over about 4000 generations (four times the mean extinction time for phage clones). …

Figure 7 with 12 supplements
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 …

Figure 7—figure supplement 1
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 …

Figure 7—figure supplement 2
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 …

Figure 7—figure supplement 3
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 …

Figure 7—figure supplement 4
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 …

Figure 7—figure supplement 5
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 …

Figure 7—figure supplement 6
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
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 …

Figure 7—figure supplement 8
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 …

Figure 7—figure supplement 9
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 …

Figure 7—figure supplement 10
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 …

Figure 7—figure supplement 11
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 …

Figure 7—figure supplement 12
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 …

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

Total population sizes equilibrate very quickly, but the total number of clones can take longer at large population sizes (high C0). The time constants inset are a measure of how quickly we expect …

Mean number of bacterial clones after t=tss 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=tss to t=5tss bacterial generations. Error bars are the standard deviation across three or more independent simulations.

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 …

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.

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 C0=104, μ=10-6, η=0.001, and e=0.95. Early time dynamics differ …

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 C0=104, μ=10-6, η=0.001, and e=0.95.

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 C0=104, μ=10-6, η=0.001, and e=0.95.

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

The first 10 trajectories that surpass nVi=1000 are shown. The simulation parameters are C0=104, μ=10-6, η=0.001, and e=0.95.

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.

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 Pnlarge in Equation 28 smoothed with a running average of window …

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 …

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

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 …

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 …

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.

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

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 es0t of Equation 31.

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

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 Equ…

Effective e=enBsnVinVinBi vs e/m across all simulations where m1 on average.

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

Four simulations with C0=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 …

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 …

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 …

nV, nB, and ν vs. η for e=0.05, approximating ν-d/c.
nV, nB, and ν vs η for e=0.15, approximating ν-c+c2-4bd2b.
nV, nB, and ν vs. η for e=0.8, approximating ν-b+b2-4ac2a.
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.

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 …

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.

Approximate phage mutation rate for e=0 vs theoretical phage mutation rate.
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 …

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 …

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 …

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 C0=104, η=10-2, μ=10-6, e=0.95, and mean m=14.5.

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 …

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

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

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 m-1, a purely phenomenological choice. Error bars are the standard deviation across three or more independent simulations.

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 m-1 and the solid line is m-1. Error bars are the standard deviation across three or more independent simulations.

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

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

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. C0=104, η=10-3, μ=10-5, e=0.95, initial m=1.

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. C0=104, η=10-3, μ=10-5, e=0.95, initial m=10.

Clan number and size over time in a simulation with C0=104, η=10-3, μ=10-5, e=0.95, initial m=1.
Clan number and size over time in a simulation with C0=104, η=10-3, μ=10-5, e=0.95, initial m=10.
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 (1.3) and m/2. Error bars are the standard deviation across three or more independent simulations.

Average clan number vs. m for all simulations that begin with one clone with μ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.

Phage infection success probability pV as a function of mutational distance between spacers and protospacers for different definitions and degrees of cross-reactivity.

Definitions are plotted for pV=0.02, e=0.95.

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.

Mean phage clone size (top) and mean bacteria clone size (bottom) relative to time of phage mutation for different definitions of pV.

These simulations begin with one initial phage clone and parameters C0=104, μ=10-6, η=0.0001, e=0.95.

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 C0=104, μ=10-6, η=0.0001, e=0.95.

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 C0=104, μ=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.

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 C0=104, μ=10-6, η=10-4, and e=0.95. Later times in this simulation are shown in Figure 55.

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 C0=104, μ=10-6, η=10-4, and e=0.95. Earlier times in this simulation are shown in Figure 54.

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.

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

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

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 (C0=104, μ=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 …

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 …

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 …

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: C0=104,η=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 …

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 …

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.

Time-shifted average immunity for three regimes of the simulation shown in Figure 54 (C0=104, μ=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 …

Time-shifted average immunity for the corresponding simulation to Figure 64 without cross-reactivity (C0=104, μ=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.
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 C0=104, μ=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 …

Average clan size across simulations with different parameters for different degrees of cross-reactivity with minit=10.

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

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 …

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 …

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.

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.

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.

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 …

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.

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.

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 …

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.

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 …

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

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.

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 Guerrero et al., 2021a.

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 …

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.

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 …

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

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

Appendix 2—figure 1
Diversity vs. η (left), μ (centre), and e (right) for C0=104.

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

Appendix 2—figure 2
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 C0=104, e=0.95, and μ=10-6 (blue). Predicted diversity in the non-coevolving model …

Appendix 2—figure 3
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 k13, independent of mutation rate (solid line). (B

Appendix 3—figure 1
Phage clone initial growth rate vs. total bacteria normalized by the initial nutrient concentration C0.

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 …

Appendix 3—figure 2
Clone size histograms (left) and cumulative distributions (right) for four different values of the spacer acquisition probability η.

In all simulations, C0=104, e=0.95, and μ=3×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 …

Appendix 3—figure 3
Bacteria and phage clone trajectories aligned to the time at which bacteria trajectories go extinct.

Bacteria trajectories are included if they reach size nBi* given by Equation 25 and all corresponding phage trajectories are plotted. In this simulation C0=3×104, e=0.8, η=10-3, and μ=10-6.

Appendix 3—figure 4
Mean time to extinction for bacterial clones after reaching size nBi* as a function of η for C0=104, e=0.95, and μ=3×10-6.

The solid line is given by Equation 139 with n=nBi*, and the dashed line is given by numerically solving Equation 137.

Appendix 3—figure 5
Predicted bacteria clone time to extinction with drift.

Measured vs. predicted mean time to extinction for bacterial clones after reaching size nBi. The predicted time to extinction is the solution with drift, given by numerically solving Equation 167.

Appendix 3—figure 6
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 …

Appendix 3—figure 7
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 nBm, where m is measured from simulations and nB is the total bacteria …

Appendix 3—figure 8
Phage clone extinction times and theoretical predictions for a simulation with parameters C0=104,η=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 nVi*=16170 as given …

Appendix 3—figure 9
Large clone extinction times from a simulation with parameters C0=104,η=0.001,e=0.95,μ=10-5.

Trajectories are counted as large if the phage clone size passes nVi*=16170, the theoretical deterministic mean phage clone size for these parameters given by Equation 26. Once a trajectory reaches nVi*, we …

Appendix 3—figure 10
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 …

Appendix 3—figure 11
Phage extinction time as a function of mean clone size.

Large phage clone mean time to extinction as a function of nVm, where m is measured from simulations and nV is the total phage population size. The solid line is given by Equation 172 without the 1+lnm

Tables

Table 1
Model reactions.
b0,i+Cg2b0,iBacterium divides
b0,iFBacterium flows out
VjFPhage flows out
FC0CNutrients flow in
CFNutrients flow out
n=0B(b0+VjαpVPn(B-n)Vj+k=1nVm+k)Interaction, phage wins,
Pn probability of n mutuant phages
b0+Vjα(1pV)(1η)b0Interaction, bacterium survives
b0+Vjα(1pV)ηbiInteraction, bacterium survives and acquires a spacer
n=0B(bi+VjαpV(i,j)Pn(Bn)Vj+k=1nVm+k)Interaction, phage wins,
Pn probability of n mutant phages
bi+Vjα(1pV(i,j))biInteraction, bacterium survives
birb0Bacterium loses spacer
Table 2
Simulation length.
C0Simulation length (bacterial generations)Steady-state start time (tss)
30010,0002000
100010,0002000
300010,0002000
10,00010,0002000
30,00015,0003000
100,00020,0004000
300,00025,0005000
1,000,00030,0006000
Table 3
Model parameters.
ParameterDescriptionValue
1gC0Bacterial doubling time41.7 min
C0Inflow nutrient concentration in102 to 10-6
Units of bacterial cell density
αPhage adsorption rate2×10-2/C0
BPhage burst size170
FChemostat flow rate0.3gC0
pVProbability of phage success0.02
for bacteria without spacers
eSpacer effectiveness0.1 to 0.95
rRate of spacer loss0.04gC0
ηProbability of spacer acquisition10-5 to 10-2
μPhage mutation rate per base per generation10-8 to 10-4
LPhage protospacer length in nucleotides30
  1. Parameter values are as above unless otherwise indicated. Representative values estimated for Streptococcus thermophilus bacteria in lab conditions.

Table 4
ν approximations.
ηeeffective
≤0.010.01 to 0.050.05 to 0.10.1 to 0.50.5 to 1
10-5-dc-dc-dc-b+b2-4ac2a-b+b2-4ac2a
10-5 to10-4-dc-dc-c+c2-4bd2b-b+b2-4ac2a
10-4 to 10-2-dc-dc-c+c2-4bd2b
10-2-dc-c+c2-4bd2b-c+c2-4bd2b
  1. Note: the drop-a solution also works for below e=0.1 for all values of η (since 2bdc20), but the -d/c solution is simpler and so preferred for the very low e range.

Additional files

Download links