1. Ecology
  2. Evolutionary Biology
Download icon

CRISPR-based herd immunity can limit phage epidemics in bacterial populations

  1. Pavel Payne  Is a corresponding author
  2. Lukas Geyrhofer
  3. Nicholas H Barton
  4. Jonathan P Bollback  Is a corresponding author
  1. University of Liverpool, United Kingdom
  2. Institute of Science and Technology Austria, Austria
  3. Technion - Israel Institute of Technology, Israel
Research Article
  • Cited 13
  • Views 2,827
  • Annotations
Cite this article as: eLife 2018;7:e32035 doi: 10.7554/eLife.32035


Herd immunity, a process in which resistant individuals limit the spread of a pathogen among susceptible hosts has been extensively studied in eukaryotes. Even though bacteria have evolved multiple immune systems against their phage pathogens, herd immunity in bacteria remains unexplored. Here we experimentally demonstrate that herd immunity arises during phage epidemics in structured and unstructured Escherichia coli populations consisting of differing frequencies of susceptible and resistant cells harboring CRISPR immunity. In addition, we develop a mathematical model that quantifies how herd immunity is affected by spatial population structure, bacterial growth rate, and phage replication rate. Using our model we infer a general epidemiological rule describing the relative speed of an epidemic in partially resistant spatially structured populations. Our experimental and theoretical findings indicate that herd immunity may be important in bacterial communities, allowing for stable coexistence of bacteria and their phages and the maintenance of polymorphism in bacterial immunity.


eLife digest

When a disease spreads through a population, it encounters certain individuals it cannot infect. If there are enough of these individuals, the epidemic stops. This phenomenon is known as ‘herd immunity’, and it occurs in many animals – for example, it plays an important role in human vaccination schemes.

While bacteria can cause disease, they are themselves targeted by viruses called ‘phages’. Bacteria can overcome this threat, and they resist phage attacks in ways that are well understood at the molecular level. However, little is known about the impact of this resistance at the scale of the population. Can herd immunity occur in bacteria? If so, what factors influence the threshold at which it will occur? In other words, what affects the minimum percentage of immune bacteria needed to stop the spread of a phage infection?

To answer these questions, Payne et al. used both experimental and mathematical methods. For the experiments, a phage and two strains of bacteria were used, one immune to the virus and one not. The two strains were combined to form several populations with different percentages of resistant bacteria, and the phage was added. How the virus could spread in these different populations was recorded. This confirmed that herd immunity does occur in bacteria and showed how the resistant bacteria influence the speed which an epidemic spreads in a population.

Building on the experiments, Payne et al. then produced a mathematical model to explore how different factors affect herd immunity. For example, the model showed that the thresholds for herd immunity can be predicted from how quickly bacteria and phages replicate. The thresholds are lower when bacteria reproduce more quickly, but higher when it is the phages that multiply faster.

The model also helps infer a formula that informs on how diseases spread in any species, such as humans. In particular, it becomes possible to predict herd immunity thresholds based on how quickly an epidemic spreads in a population where few people are vaccinated. Future research is needed to adapt the formula to the specific factors that shape disease outbreaks in humans. Ultimately, this could help policymakers design strategies to deal with infectious diseases, such as yearly outbreaks of the flu.



The term ‘herd immunity’ has been used in a variety of ways by different authors (see Fine et al., 2011). Here, we define it as a phenomenon where a fraction of resistant individuals in a population reduces the probability of transmission of a pathogen among the susceptible individuals. Furthermore, if the fraction of resistant individuals in a population is sufficiently large the spread of a pathogen is suppressed. Experimental research into the phenomenon has focused mostly on mammals (Jeltsch et al., 1997; Mariner et al., 2012), birds (van Boven et al., 2008; Meister et al., 2008), and invertebrates (Konrad et al., 2012; Wang et al., 2013). In human populations the principles of herd immunity were employed to limit epidemics of pathogens through vaccination programs (Fine et al., 2011), which in the case of smallpox lead to its eradication between 1959 and 1977 (Fenner, 1993).

Alongside advances in vaccination programs, the formalization of a general theory of herd immunity was developed. The theory is based on a central parameter, R0, which describes the fitness of the pathogen, as measured by the number of subsequent cases that arise from one infected individual in a population (for a historical review of R0 see [Heesterbeek, 2002]). Thus, R0 indicates the epidemic spreading potential in a population. Given R0 the herd immunity threshold is defined as,

(1) H=R0-1R0,

which determines the required minimum fraction of resistant individuals needed to halt the spread of an epidemic. R0 and subsequently also H are affected by the specific details of transmission and the contact rate among individuals (Grassly and Fraser, 2008). Many theoretical studies have addressed the influence of some of these details, in particular maternal immunity (Anderson and May, 1992-08), age at vaccination (Anderson and May, 1982; Nokes and Anderson, 1988), age related or seasonal differences in contact rates (Schenzle, 1984; Anderson and May, 1985; Yorke et al., 1979), social structure (Fox et al., 1971), geographic heterogeneity (Anderson and May, 1984; Lloyd and May, 1996; Real and Biek, 2007), and the underlying contact network of individuals (Ferrari et al., 2006).

Interestingly, little work has focused on the potential role of herd immunity in microbial systems which contain a number of immune defense systems and have an abundance of phage pathogens. These defenses vary in their potential to provide herd immunity as they target various stages of the phage life cycle, from adsorption to replication and lysis. Early defense mechanisms include the prevention of phage adsorption by blocking of phage receptors (Nordström and Forsgren, 1974), production of an extracellular matrix (Hammad, 1998; Sutherland et al., 2004), or the excretion of competitive inhibitors (Destoumieux-Garzón et al., 2005). Alongside these bacteria have evolved innate immune systems that target phage genomes for destruction. These include host restriction-modification systems (RMS) (Blumenthal and Cheng, 2002), argonaute-based RNAi-like systems (Swarts et al., 2014), and bacteriophage-exclusion (BREX) systems (Goldfarb et al., 2015). In addition to innate systems, bacteria have evolved an adaptive immune system called CRISPR-Cas (clustered regularly interspaced short palindromic repeat) (Sorek et al., 2013). In order for any of these immune systems to provide herd immunity, they must prevent further spread of the pathogen. Therefore, unless the phage particles degrade in the environment at a timescale comparable to the phage adsorption rate, the immune system must provide a ‘sink’ for the infectious particles reducing the average number of successful additional infections below one. Unlike the early defense mechanisms that may simply prevent an infection but not the further reproduction of infectious particles, the RMS, BREX, argonaute-based RNAi-like, and the CRISPR-Cas systems degrade foreign phage DNA after it is injected into the cell, and thus continue to remove phage particles from the environment, which increases their potential to provide herd immunity. In order for herd immunity to arise, the population must also be polymorphic for immunity, which can be achieved if immunity is plasmid borne. In addition to this, the CRISPR-Cas system is unique in that it is adaptive allowing cells to acquire immunity upon infection (see Figure 1A,B and C), which can lead to polymorphism in immunity even if the system is chromosomal.

Mechanism of CRISPR/Cas type II immunity.

The CRISPR/Cas system provides immunity to phages and its main features can be described by three distinct stages. (A) Acquisition. When a cell gets infected by a phage, a protospacer on the invading phage DNA (indicated as a red bar) is recognized by Cas1 and Cas2. The protospacer is cleaved out and ligated to the leader end (proximal to the Cas genes) of the CRISPR array as a newly acquired spacer (red diamond). (B) Processing. The CRISPR array is transcribed as a Pre-crRNA and processed by Cas9 (assisted by RNaseIII and trans–activating RNA, not shown) into mature crRNAs. (C) Interference. Mature crRNAs associate with Cas9 proteins to form interference complexes which are guided by sequence complementarity between the crRNAs and protospacers to cleave invading DNA of phages whose protospacers have been previously incorporated into the CRISPR array. (D) A truncated version of the CRISPR system on a low copy plasmid, which was used in this study lacks cas1 and cas2 genes and was engineered to target a protospacer on the T7 phage chromosome to provide Escherichia coli cells with immunity to the phage. The susceptible strain contains the same plasmid except the spacer does not target the T7 phage chromosome.


In addition to immune system-specific factors, the reproductive rate of phage depends strongly on the physiology of the host bacterium (Hadas et al., 1997), and the underlying effective contact network which may vary greatly in bacterial populations depending on the details of their habitat. Thus, herd immunity will be influenced by the physiological state of the bacteria and the mobility of the phage in the environment through passive diffusion and movement of infected individuals. Taken together these details call into question the applicability of the traditional models of herd immunity from vertebrates to phage-bacterial systems. Thus, experimental investigation and further development of extended models that take into account the specifics of microbial systems are required.

To investigate under which conditions herd immunity may arise in bacterial populations, we constructed an experimental system consisting of T7 phage and bacterial strains susceptible and resistant to it. Our experimental system can be characterized by the following features. First, we used two strains of Escherichia coli, one with an engineered CRISPR-based immunity to the T7 phage, and the other lacking it (Figure 1D). Second, we examined the dynamics of the phage spread in different environments – spatially structured and without structure. Furthermore, we developed and analyzed a spatially explicit model of our experimental system to determine the biologically relevant parameters necessary for bacterial populations to exhibit herd immunity.


Properties of resistant individuals

We engineered a resistant E. coli strain by introducing the CRISPR-Cas Type II system from Streptococcus pyogenes with a spacer targeting the T7 phage genome (see Material and Methods). We further characterized the ability of the system to confer resistance to the phage. We find a significant level of resistance as measured by the probability of cell burst when exposed to T7 (Figure 2A). However, resistance is not fully penetrant as approximately 1 in 1000 resistant cells succumb to infection. In addition, we observe that as phage load increases (multiplicity of infection, MOI) the probability that a cell bursts increases (Figure 2A). In order to determine the herd immunity threshold in our experimental system, we constructed the resistant strain such that upon infection the cell growth is halted, yet the cell still adsorbs and degrades phages (Figure 2B,C). This feature is important as it prevents the action of frequency dependent selection which in naturally growing populations will favor the resistant strain until its frequency reaches the herd immunity threshold. Thus, in our system if the frequency of the resistant strain is below the herd immunity threshold, the resistant cells remain below the threshold and are unable to stop the epidemic and the whole population collapses. In contrast, if the frequency of resistant individuals in the population is above the herd immunity threshold, the resistant individuals provide complete herd immunity and the population survives. These properties allow us to quantify the expanding epidemic in both liquid media and on bacterial lawns (without and with spatial structure, respectively) using high throughput techniques. Specifically, it allows us to control for the complex dynamics of the system arising from frequency dependent selection and simultaneous changes in the physiological states of the cells (growth rates depending on the nutrient concentrations) and phage (burst size, latent period depending on the cell’s physiology).

Efficiency of bacterial resistance.

(A) The probability that a resistant cell bursts, relative to a susceptible cell, at three different initial multiplicities of infection (MOI). The probability that a resistant cell bursts at MOI 1000 is significantly higher than at MOI 10 (p = 0.019, t4 = 3.031) or at MOI 100 (p = 0.022, t5 = 2.674). The error bars show the standard deviations from the mean. Note that this measure is not a widely used ’efficiency of plating’ but it determines the probability of burst of single resistant cells (see Materials and methods for details). (B) The number of colony forming units (CFUs) post phage challenge (see Materials and methods). The mean number of CFUs after the bacterial cultures were exposed to the phage is not significantly different between susceptible and resistant strains at MOI 10 (p = 0.239, t22 = 0.721) and (C) at MOI 100 (p = 0.27, t30 = 1.124), indicating that the resistant cells’ growth is halted after the cells are infected by a phage. The error bars show the standard deviations from the mean. There were no detectable CFUs in either susceptible or resistant cell cultures at MOI 1000. It should be noted that the indicated MOI values do not correspond to the average number of phages that adsorb to cells in the experiments. For MOI 10 we estimated the mean number of phages per cell as 0.229 and for MOI 100 as 0.988 (see Materials and methods for details). It was impossible to determine the mean for MOI 1000 as there were no detectable CFUs under such conditions. The data presented in this figure can be found in Figure 2—source data 1.

Figure 2—source data 1

Efficiency of bacterial resistance.


It should be noted that our model does not reflect this artificial property – it assumes that resistant bacteria keep growing after successfully overcoming a phage infection (see Equation (2d)). This discrepancy, however, does not affect the model prediction of the herd immunity threshold in our experimental system for the following reason: time scale of an epidemic spread through a population (double exponential phage growth) is substantially shorter than the time scale of bacterial population growth (exponential growth). Therefore, whether or not an epidemic is established does not depend on later dynamics of frequencies of resistant and susceptible individuals in the population, it only depends on the initial conditions. Similarly, the model correctly captures the dynamics of an epidemic in spatially structured populations as the phage spreads radially and in every time-point the epidemic front encounters a naive population with a constant ratio of resistant to susceptible individuals.

Herd immunity in populations without spatial structure

To understand the influence of spatial population structure, or lack thereof, we first measured the probability of population survival (i.e., whether the cultures are cleared or not) in well mixed liquid environments (no spatial structure) consisting of differing proportions of resistant to susceptible individuals and T7 phage. When the percentage of resistant individuals is in excess of 99.6% all 16 replicate populations survive a phage epidemic (i.e., show no detectable difference in growth profiles to the phage free controls; Figure 3). Populations with 99.2% and 98.4% resistant individuals show intermediate probabilities of survival – 10 out of 16 replicate populations and 4 out of 16 replicate populations survive, respectively (Figure 3). The likely explanation as to why some populations survive and others collapse is due to the stochastic nature of phage adsorption after inoculation: If the population composition is close to the herd immunity threshold a stochastic excess of phage particles adsorbing to susceptible cells may trigger an epidemic, whereas if chance increases the number of phages adsorbing to resistant individuals, the epidemic is suppressed. However, when populations have fewer than 96.9% resistant individuals all 16 replicate populations fail to survive and collapse under the epidemic (Figure 3).

Fraction of surviving populations at 18h post phage infection.

Bacterial populations consisting of various fractions of resistant to susceptible individuals infected with 50 phages, corresponding to a multiplicity of infection (MOI) of 10-4, designed to resemble an epidemic initiated by the burst size from one infected individual (see Table 2 for burst size estimates). Each population phage challenge is replicated 16 times. The solid dark green line shows the model prediction, Equation 4, for the herd immunity threshold (H), given latent period (λ), bacterial growth rate (α), and phage burst size (β). Shaded area indicates ±1 standard deviation. The data presented in this figure can be found in Figure 3—source data 1.

Figure 3—source data 1

Fraction of surviving populations at 18 hr post phage infection.


As mentioned in the introduction, phage and bacterial physiology may affect the herd immunity threshold. To test this we altered bacterial growth by reducing the concentration of nutrients in the medium by mixing LB broth with 1X M9 salts in different ratios (Figure 4), which concurrently alters the T7 phage’s latent period and burst size (Figure 5A,B and Table 1). Indeed, we observe as bacterial growth rates decline the fraction of resistant individuals necessary for population survival decreases (Figure 5C). When the populations are grown in a 50% diluted growth medium, the fraction of resistant individuals required for a 100% probability of survival is 99.2%; when the fraction of resistant individuals is 75% or less populations go extinct. In a 20% growth medium the fraction of resistant individuals required for survival decreases to 96.9%, while the fraction when all replicates collapse to 50%.

Measuring bacterial growth without phage.

(A) Trajectory of population size on agar plates over time. For modeling, we assume two states of growth (dashed brown curve): first, the bacterial population grows exponential until the time Tdepl, when nutrients are depleted. From this time on, growth rate is assumed to be zero and the population saturates at a maximal size Bfinal. Experimental observations fit this proposed growth curve to a very good extent. After all, half of all nutrients are used up in the last generation indicating that the switch between growth and no-growth should be fast. (B) Growth rates of bacteria in diluted medium follow closely Monod’s empirical law, given by expression Equation 9. Fit parameters are found to be αmax0.720h-1 and Kc0.257 (with the latter in dimensionless units as dilution of LB medium), see also Table 1. The data presented in this figure can be found in Figure 4—source data 1.

Figure 4—source data 1

Bacterial growth on soft agar plates (tab Figure 4A) and bacterial growth in LB medium of various dilutions (tab Figure 4B).

Herd immunity threshold in liquid culture as a function of bacterial growth.

(A) Phage burst size (β) change as a function of nutrient concentrations. (B) Latent period (λ) increase across the range of nutrient concentrations. Values for β and λ are given in Table 2. (C) Population survival analysis upon phage challenge as a function of the fraction of resistant cells and the intrinsic growth rate (nutrient availability, N). Bacteria survive the phage infection (full circles), collapse (empty circles), or exhibit both outcomes (circled dots) in the 16 to 18 replicates, done in 3 independent batches. Light green errorbars at investigated dilutions of LB show the expected value and its standard deviation of H(α), Equation 5, with standard error propagation of the measured β, λ and α. In order to interpolate herd immunity to dilutions N not probed in experiments (dark green line), we use a second order polynomial in N to fit the data for both β/λ and β, which excellently matches the average measurements (a naive linear fit displays non-negligible deviations and non-sensical negative values). In addition, the dependence α=α(N) is obtained by numerically inverting the Monod growth rate dependence, see Equation 9. The data presented in this figure can be found in Figure 5—source data 1 and Figure 5—source data 2.

Figure 5—source data 1

Figure 5A B source data: Phage burst sizes and latent period in different dilutions of the growth medium.

Figure 5—source data 2

Figure 5C: Fraction of surviving populations in different dilutions of the growth medium.

Table 1
Estimated parameters for bacterial growth using Monod kinetics.

Undiluted LB medium (N=1) is assumed to have 15mg/ml nutrients (10mg/ml Tryptone, 5mg/ml yeast extract). The full dataset is shown in Figure 4.

αmax0.720(± 0.011)[h-1]
Kc0.257(± 0.012)Dilution N of LB [01]

From the experimental observations of the herd immunity threshold values we infer the phage R0 using Equation 1. In an undiluted growth medium the phage R0 falls between 32 and 256 and decreases to between 4 and 128 in 50% and between 2 and 32 in 20% nutrient medium. These data indicate that bacterial populations can exhibit herd immunity in homogeneous liquid environments. However, bacteria typically live in spatially structured environments such as surfaces, biofilms or micro-colonies, therefore we extended our experiments to consider the potential impact of spatially structured populations.

Herd immunity in spatially structured populations

In order to discern the role, if any, spatial structure plays in herd immunity we conducted a set of experiments in spatially structured bacterial lawns on agar plates. Spatially structured bacterial populations provide a more fine grained measure of herd immunity, compared to the population survival assays done in liquid culture. On bacterial lawns, phages spread radially from a single infectious phage particle and the radius of plaque growth on different proportions of resistant to susceptible individuals can be easily quantified. In addition, these data allow for estimating the speed of the epidemic wave front in these different regimes using real-time imaging (Figure 6A).

Properties of expanding phage epidemics on bacterial lawns.

(A) Example of plaque morphology and size change over 48 hours for populations with 50% resistant cells (top) and a control with 100% susceptible cells (bottom). (B) Mean plaque size area through time. Colors indicate the different fraction of resistant individuals (color coding as in panel C). Note the distinct two phases of plaque growth – initially, phage grow fast with exponentially growing bacteria but slow once the nutrients are depleted (10 hr). The plaque radius is reduced, relative to 100% susceptible population, even when only a small fraction of resistant individuals are in the population. (C) Final plaque radius at 48 hpi. Green line shows the prediction from the model for the plaque radius r. Grey numbers indicate the number of plaques measured. Error bars indicate the standard deviations. The data presented in this figure can be found in Figure 6—source data 1.

Figure 6—source data 1

Plaque radii for all population compositions and time points.


We observe a decline in the number of plaque forming units (see Appendix 2—figure 1) and a significant decrease in final plaque sizes as the proportion of resistant individuals in the populations increases (Figure 6B,C). A reduction in the final plaque size compared to a fully susceptible population was statistically significant with as few as 10% resistant individuals in a population (p=0.004, t53 = 2.744). In order to determine the effect of resistant individuals during the earlier phase of bacterial growth (until the bacterial density on the agar plate reaches saturation; Figure 4A), we analyze the velocities of plaque growth between 0 and 24 hr post inoculation (hpi). We find that the speed is significantly reduced after 11hpi when the population consists of as few as 10% of resistant individuals (p=0.0317, t32 = 1.923). As the fraction of resistant individuals further increases, the speed declines significantly at earlier and earlier time points: 6hpi with 20% (p=0.0392, t62 = 1.79), and 5.67hpi with 30% (p=0.0286, t53 = 1.943). In fact, when the fraction of resistant individuals exceeds 40%, the reduction in the speed of the spread is statistically significant immediately after the plaques are visually detectable (Figure 7). It should be noted that all populations with such low percentages of resistant individuals in liquid environment collapsed, indicating that spatial structure plays a role in herd immunity.

Speed of phage epidemic expansion on bacterial lawns.

(A) Speed of expanding phage epidemics for all population compositions is initially high, before it drops once nutrients are depleted at around 10hpi (hours post infection). (B) Plaque speed significance. Comparing speeds of plaque spread with the 100% susceptible control. Linear regression of a sliding window spanning 4 hours of the radius sizes was calculated for all individual plaques and all compositions of the populations between t0 and t24. Slopes of the linear regressions for all compositions of the populations were compared using a two-sided heteroscedastic t-test against the 100% susceptible dataset. The data presented in this figure can be found in Figure 7—source data 1.

Figure 7—source data 1

Speed of plaque expansion in populations consisting of varying proportions of resistant to susceptible bacteria.


The results presented in this and the previous section would allow us to use Equation 1 to infer a value for R0 from the observed threshold between surviving and collapsing bacterial populations, Figures 3,5. We also observe that herd immunity is strongly influenced by spatial organization of the population, Figure 6. How the exact value of H (and subsequently the ‘classical’ epidemiological parameter R0) is affected by various factors such as bacterial growth rate, phage burst size and latent period is, however, difficult to resolve experimentally. Similarly, quantification of the effect of spatial structure is hardly achievable solely by experimental investigation. In order to disentangle the roles of all the factors mentioned above, we proceed with development and analysis of a mathematical model of the experimental system.

Modelling bacterial herd immunity

We developed a model of phage growth that takes several physiological processes into account: bacterial growth during the experiment, bacterial mortality due to phage infection, and phage mortality due to bacterial immunity. Furthermore, we use the previously reported observation that phage burst size β and latent period λ depend strongly on the bacterial growth rate α (see Table 1).

The main processes in our model system can be defined by the following set of reactions,


Susceptible (Bs) and resistant (Br) cells grow at a rate α (no significant difference in growth rate between strains, α(Bs)=0.551±0.045h-1, α(Br)=0.535±0.023h-1,t70=1.867,p=0.066), Equation 2, by using an amount y of the nutrients N. Phage infection first involves adsorption to host cells, Equation 2c and Equation 2d, with the adsorption term A specified below. Infected susceptible bacteria produce on average β phage with a rate inversely proportional to the average latency λ. In contrast, resistant bacteria either survive by restricting phage DNA via their CRISPR-Cas immune system or – less likely – succumb to the phage infection. However, when the MOI is large even resistant cells become susceptible to lysis resulting in the release of phage progeny (see Figure 2) (Westra et al., 2015; Chabas et al., 2016).

In our system, bacteria eventually deplete the available nutrients, N(t>Tdepl)=0, resulting in the cessation of growth. This decline in bacterial growth affects phage growth – latency increases and burst size decreases, such that phage reproduction declines dramatically (see Table 2). We define the critical time point at which cells transition from exponential growth to stationary phase as,

(3) Tdepl1αlog(BB0).
Table 2
Estimated parameters for phage growth.

See also Figure 5A,B.

MediumDilutionLatent periodBurst sizeBurst size/hour
LB 00.0101.1 (±10.9)3.0 (±1.9)1.8 (±1.1)
LB 200.243.4 (±3.9)3.0 (±1.9)16.6 (±6.0)
LB 500.540.0 (±3.0)35.6 (±16.4)53.4 (±24.9)
LB 1001.036.1 (±6.1)85.6 (±47.3)142.1 (±82.1)

Here, B0 and B are the initial and final bacterial densities, respectively. In the initial exponential growth phase, our estimates from experimental data for growth parameters are α=0.63h-1, β=85.6phages/cell and λ=0.60h, for bacteria and phages, respectively (Tables 1 and 2). After time Tdepl, bacterial growth rate is set to zero (α=0) and phage growth is reduced to βdepl=3.0phages/cell and λdepl=1.69h. Such a two state model – constant growth rate while nutrients are present and no growth after depletion – describes the observed population trajectories in experiments sufficiently well (see Figure 4).

Modelling herd immunity in populations without structure

An important parameter for estimating herd immunity is the fraction S of susceptible bacteria in the population. As a first estimate, a phage infection spreads in well mixed bacterial cultures if βS>1, which leads to a continuous chain of infections: the product of burst size β of phage particles and the probability S of infecting a susceptible host has to be larger than one. As a first approximation, one could identify R0 with the burst size β, which is compatible with the observed herd immunity thresholds when inverting Equation 1.

However, the growing bacterial population could outgrow the phage population if the former reproduces faster (e.g., in the case of RNA coliphages, van Duin, 1988), which introduces deviations from the simple relation between R0 and H as shown in Equation 1. We capture this dynamical effect in a correction to the previous estimate as βS>1+λα (see Materials and methods): more phages have to be produced for the chain of infections to persist in growing populations. The correction λ1/α is the ratio of generation times of phages over bacteria – usually, such a correction is very small for non-microbial hosts and can be neglected. Ultimately, herd immunity is achieved if the threshold defined by H=1-Sc is exceeded, with Sc the critical value in the inequality above. Rearranging, we obtain an expression for the herd immunity threshold

(4) H=β-1-λαβ.

This estimate of H coincides to a very good extent with the population compositions of susceptible and resistant bacteria where we observe the transition from surviving and collapsed populations in experiments (see Figure 3). Moreover, simulations presented in the Appendix (section Simulation of recovery rate) show a range in the bacterial population composition with non-monotonic trajectories for Bs and Br (see Appendix 1—figure 1B), which is comparable to the range in composition we find in both outcomes, that is, some surviving and some collapsing populations in experiments. For such parameter choices, stochastic effects could then decide the observed fates of bacteria. As presented above, the herd immunity threshold changes when the bacterial cultures grow in a diluted growth medium. In a set of independent experiments we measured bacterial growth rate α, phage burst size β and phage latent period λ under such conditions (see Figure 4B and Table 2). From these data we estimated the dependence of the phage burst size on the bacterial growth rate, β(α), using a numerical quadratic fit (Figure 5A). Similarly, we estimated the dependence of the phage latent period on the bacterial growth rate, λ(α) (Figure 5B). Using these estimates we calculated the expected growth rate–dependent herd immunity threshold

(5) H(α)=β(α)-1-λ(α)αβ(α),

which gives a very good prediction of the shift in the herd immunity threshold to lower values for slower growing populations (green line in Figure 5C). This verification of our model shows that it correctly captures the dependence of the herd immunity threshold on bacterial and phage growth parameters.

The deviations from the herd immunity threshold depicted by the green area in Figure 3 and green error bars in Figure 5C are derived from uncertainty in measurements in β, λ and α. The inherent stochasticity of the adsorption process thus provides additional uncertainty, which is not captured by the depicted error bars. This additional stochasticity can explain wider transition zone in experiments with slower growing populations (dilution 0.5 and 0.2), because the fate of the population is more prone to stochastic effects as the phage replication rate is slower than in a fast growing population. This stochastic effect might be reduced by larger phage inocula. This could, however, also shift the observed transition between collapsing and surviving populations towards higher frequencies of resistant bacteria (and away from the actual herd immunity threshold) as protection by the immune system is less effective with increasing number of phages per cell (see Figure 2A).

Modelling herd immunity in spatially structured populations

The dynamics of phage spread differ between growth in unstructured (e.g., liquid) and structured (e.g., plates) populations. In order to quantify the effect of spatial structure in a population, we extend our model to include a spatial dimension. In structured populations growth is a radial expansion of phages defined by the plaque radius r and the expansion speed v, for which several authors have previously derived predictions (Kaplan et al., 1981; Yin and McCaskill, 1992; You and Yin, 1999; Fort and Méndez, 2002a; Ortega-Cejas et al., 2004; Abedon and Culler, 2007; Mitarai et al., 2016).

We assume phage movement can be captured by a diffusion process characterized with a diffusion constant D, which we estimate in independent experiments as D=1.17(±0.26)102 mm2/h (see Materials and methods, Figure 8). However, we assume that only phages disperse and bacteria are immobile as the rate of bacterial diffusion does not influence the expanding plaque on timescales relevant in the experiment. Adsorption of phages on bacteria is modeled with an adsorption constant δ.

Estimating diffusion constant of phages.

(A), (B) Phage are slowly expanding on agar which can be observed via their fluorescence. Pictures are taken 5h apart. (C) The diffusion constant D can be estimated as best-fit parameter in a heat kernelK(D) propagates the fluorescence profile L(t) at time t forward (via a convolution to “smear” out the signal) to the profile L(t+Δt) at the next measured time point. The difference between the expected change and the actual profile is quantified as total squared deviation, see Equation 10, which we minimize to obtain D. Consequently, we can estimate the diffusion constant as D1.1710-2mm2/h. The green line uses this estimated parameter D and shows the change between the profile at t=10h (orange line) and the profile at t=15h (light brown line), assuming diffusive spread of phages. See Materials and methods for more information.


Taking these considerations together, allows to write a reaction-diffusion dynamics for growth of phages P on the growing bacterial population as

(6) tP=Dx2P+δ(βS-1-λα)P.

The first term accounts for the diffusive spread of phages, while the second term describes phage growth. This second term includes the correction λα which arises due reproduction of bacteria, derived in the unstructured liquid case.

The spreading infection will sweep across the bacterial lawn with the following speed

(7) v=2DδβS-1-λα,

which is computed in more details in the Materials and methods. This expression Equation 7 indicates that the population composition crucially influences the spreading speed at much lower fractions of resistant bacteria than the herd immunity threshold Equation 4, where phage expansion stops completely. Consequently, the resulting plaque radius r decays with increasing fractions of resistants and reaches zero at H. A prediction for r can be obtained by integrating Equation 7 over time.

In our (simplified) model, time-dependence of the speed only enters via the fraction of susceptibles S, which is assumed to stay at the initial S0 value until it encounters the epidemic wave of phages. Furthermore, we use the experimental observation that plaque expansion ceases upon depletion of nutrients, coinciding with a cessation of bacterial growth. This leads to the approximation rvTdepl, with Tdepl given by Equation 3. Using this expression we estimated the adsorption constant δ from the growth experiments as it is difficult in practice to measure on plates. The green line in Figure 6B is the best fit, yielding the value δ=4.89(±0.19)10-2bacteria/phageh for the adsorption constant.

Our results for spatially structured populations allows us to speculate on a general epidemiological question: If an infection is not stopped by herd immunity in a partially structured population, by how much is its spread slowed down? By generalizing Equation 7 we can derive a relative expansion speed, compared to a fully susceptible population,

(8) vrel=1-1-SH.

This expression, Equation 8, is devoid of any (explicit) environmental conditions, which are not already contained in the herd immunity threshold H itself. Thus, it could apply to any pathogen-host system. Ultimately, this relative speed approaches zero with a universal exponent of 1/2, when the fraction of resistant individuals 1-S approaches the herd immunity threshold H. However, a few caveats exist when using Equation 8, as several conditions have to be fulfilled: Obviously, a pathogen is expected not to spread in a population exhibiting complete herd immunity – the relative speed should only hold for populations below the herd immunity threshold. Moreover, if dispersal cannot be described by diffusion, but rather dominated by large jumps (Hallatschek and Fisher, 2014), the diffusion approach we used for traveling waves is not applicable, and thus also renders Equation 8 inadequate.

An increase in the number of long range jumps of phages can be considered as a transition between the two cases we treated here – spatially explicit dynamics on plates and completely mixed populations in liquid culture, respectively. Potential long range jumps of phages can be mediated by host cells moving distances that the phages cannot achieve on their own. In such cases, dispersal of the phages is a convolution of movement of their hosts with their own ability to spread locally. These long range jumps would therefore increase the overall expansion speed and area of the epidemic. We expect that in our setup bacterial motility does not substantially contribute to phage spread because (i) bacteria become motile only in late exponential/early stationary phase (Amsler et al., 1993) when phage reproduction drops to very low levels, and (ii) the soft agar concentration used in our experiments (0.525%) effectively blocks bacterial motility (Croze et al., 2011). However, we would not expect that long range jumps change the herd immunity threshold H(α) itself. Spread of pathogens still stops when the fraction of susceptible hosts S is small such that βS<1+λα, and will continue as long as βS>1+λα is fulfilled.


The spread of a pathogen may be halted or slowed by resistant individuals in a population and thus provide protection to susceptible individuals. This process, known as herd immunity, has been extensively studied in a wide diversity of higher organisms (Jeltsch et al., 1997; Mariner et al., 2012; van Boven et al., 2008; Meister et al., 2008; Konrad et al., 2012; Wang et al., 2013). However, the role of this process has largely been ignored in microbial communities. To delve into this we set out to determine under what conditions, if any, herd immunity might arise during a phage epidemic in bacterial populations as it could have profound implications for the ecology of bacterial communities.

We show that herd immunity can occur in phage-bacterial communities and that it strongly depends on bacterial growth rates and spatial population structure. Average growth rates of bacteria in the wild have been estimated as substantially slower than in the laboratory (generation time is 7.4 fold longer [Gibson et al., 2017]), a condition that we have shown to facilitate herd immunity. Furthermore, bacterial populations in the wild are also highly structured, as bacteria readily form micro-colonies or biofilms (Hall-Stoodley et al., 2004) and grow in spatially heterogeneous environments such as soil or the vertebrate gut (Fierer and Jackson, 2006), a second condition we have shown to facilitate herd immunity. These suggest that herd immunity may be fairly prevalent in low nutrient communities such as soil and oligotrophic marine environments.

In an evolutionary context, herd immunity may also impact the efficacy of selection as the selective advantage of a resistance allele will diminish as the frequency of the resistant allele in a population approaches the herd immunity threshold, H. This has two important implications. First, while complete selective sweeps result in the reduction of genetic diversity at linked loci, herd immunity may lead to only partial selective sweeps in which some diversity is maintained. Second, herd immunity has a potential to generate and maintain polymorphism at immunity loci, as has been shown for genes coding for the major histocompatibility complex (MHC) (Wills and Green, 1995). Polymorphism in CRISPR spacer contents have been demonstrated in various bacterial (Tyson and Banfield, 2008; Sun et al., 2016; Kuno et al., 2014) and Archaeal (Held et al., 2010) populations and communities (Pride et al., 2011; Zhang et al., 2013; Andersson and Banfield, 2008). While these studies primarily explain polymorphisms in CRISPR spacer content as a result of rapid simultaneous independent acquisition of new spacers, we suggest that observed polymorphisms may result from frequency-dependent selection on resistance loci arising from herd immunity. In such a case, herd immunity is likely to maintain existing polymorphism in CRISPR spacer content in 1-H fraction of the population, unless the current major variant goes to fixation due to drift. However, considering the large population sizes of bacteria, drift is unlikely to have a strong effect, allowing herd immunity to maintain a large fraction of immunity polymorphism.

It has also been suggested that herd immunity might favor coexistence between hosts and their pathogens (Hamer, 1906), which can lead to cycling in pathogen incidence and proportions of resistant and susceptible individuals over time, e.g., in measles before the era of vaccination (Fine, 1993). This cycling is caused by the birth of susceptible individuals, which, once their proportion exceeds the epidemic threshold (1-H), lead to recurring epidemics. CRISPR-based immunity is, however, heritable meaning that descendants of resistant bacteria remain resistant. One might speculate that analogous cycling in phage epidemics may occur if immunity is costly. In turn, a computer simulation study of coevolution of Streptococcus thermophilus and its phage found both cycling and stable coexistence of different CRISPR spacer mutants and phage strains (Childs et al., 2014). The extent to which herd immunity facilitates maintenance of CRISPR spacer polymorphism and coexistence with phage requires further experimental and theoretical investigation.

We also developed a mathematical model and show how the herd immunity threshold H (Equation 4) depends on the phage burst size β and latent period λ, and on the bacterial growth rate α. This dependence arises as phages have to outgrow the growing bacterial population, as host and pathogen have similar generation times in our microbial setting. In addition to these parameters, we also describe how the speed v (Equation 7) of a phage epidemic in spatially structured populations depends on phage diffusion constant D, phage adsorption rate δ, and the fraction of resistant and susceptible individuals in the population. All of which are likely to vary in natural populations. We also derived the relative speed of spread for partially resistant populations, as measured relative to a fully susceptible population, and show that it can be parametrized solely with the herd immunity threshold H (Equation 8). This relative speed of the spread of an epidemic should be applicable to any spatially structured host population where the spread of the pathogen can be approximated by diffusion. Both our experiments and the modelling show that even when the fraction of resistant individuals in the population is below the herd immunity threshold the expansion of an epidemic can be substantially slowed, relative to a fully susceptible population.

In conclusion, we have presented an experimental model system and the connected theory that can be usefully applied to both microbial and non-microbial systems. Our theoretical framework can be useful for identifying critical parameters, such as H (and to some extent R0), from the relative speed of an epidemic expansion in partially resistant populations so long as the process of pathogen spread can be approximated by diffusion. This approximation has been shown to be useful in such notable cases as rabies in English foxes (Murray et al., 1986), potato late blight (Scherm, 1996), foot and mouth disease in feral pigs (Pech and McIlroy, 1990), and malaria in humans (Gaudart et al., 2010).

Materials and methods

Key resources table
Table of key strains, reagents and software used in this study.
Reagent type (species)
or resource
DesignationSource or referenceIdentifiersAdditional information
gene (Streptococcus
pyogenes SF370)
cas9National Center for
Biotechnology Information
Gene symbol SPy_1046
strain, strain background
(Escherichia coli)
E. coli K12 MG1655Own collectionNA
strain, strain background
(Bacteriophage T7)
E. coli bacteriophage T7ATCC CollectionATCC:BAA-1025-B2;
recombinant DNA reagentpCas9Addgene Vector DatabaseAddgene:42876;
pCas9 plasmid was a gift
from Luciano Marraffini
recombinant DNA reagentpCas9T7resistantthis paperNAPlasmid derived from pCas9
commercial assay or kitPureYield Plasmid
Miniprep System
chemical compound, drugChloramphenicolSigma-AldrichSigma-Aldrich:C0378-5G;
software, algorithmPerkinElmer Volocity v6.3RRID:SCR_002668Volocity 3D Image
Analysis Software
software, algorithmFiji v1.0doi: 10.1038/nmeth.2019RRID:SCR_002285Image processing
package of ImageJ
software, algorithmRStudio 1.0.153RRID:SCR_000432Software for the R
statistical computing
software, algorithmPython 3.6.3RRID:SCR_008394Python programming language
software, algorithmModel source codedoi: 10.5281/zenodo.1038582RRID:SCR_004129Zenodo repository

Experimental methods

Engineering resistance

Request a detailed protocol

Oligonucleotides AAACTTCGGGAAGCACTTGTGGAAG and AAAACTTCCACAAGTGCTTCCCGAA were ordered from Sigma-Aldrich, annealed and ligated into pCas9 plasmid (pCas9 was a gift from Luciano Marraffini, Addgene plasmid #42876) carrying a Streptococcus pyogenes truncated CRISPR type II system and conferring chloramphenicol resistance. For the detailed protocol see (Jiang et al., 2013). The oligonucleotides were chosen so that the CRISPR system targets an overlap of phage T7 genes 4A and 4B. Therefore, the CRISPR system allows the gene 0.7, coding for a protein which inhibits the RNA polymerase of the cell, to be expressed before the T7 DNA gets cleaved (García and Molineux, 1995). The subsequent growth of the cells is halted and phage replication is inhibited. The plasmid was electroporated into Escherichia coli K12 MG1655 (F- lambda- ilvG- rfb-50 rph-1). The T7 wildtype phage was used in all experiments.

Efficiency of the CRISPR-Cas system

Request a detailed protocol

Efficiency of the engineered CRISPR-Cas system was tested using the following protocol: Overnight cultures grown in LB containing 25 μg/ml chloramphenicol were diluted 1 in 10 in the same medium, cells were infected with the T7 phage (MOI 10, 100, and 1000) and incubated for 15min in 30C. Cells were spun down for 2 min in room temperature at 21130g. Supernatant was discarded and the cell pellet was resuspended in 950 μl of 1X Tris-HCl buffer containing 0.4% (227μM) ascorbic acid pre-warmed to 43C and incubated in this temperature for 3 min to deactivate free phage particles (Murata and Kitagawa, 1973). Cultures were serially diluted and plated using standard plaque assay protocol on a bacterial lawn of susceptible cells to detect bursting infected cells. The supernatant was tested for free phage particles, which were not detected in the corresponding dilutions used for plaque counting. Each experiment was replicated three or four times (MOI 10 three times, MOI 100 four times and MOI 1000 three times) while samplings from each treatment were performed in quadruplicates. The probability that a resistant cell bursts was calculated as a ratio of bursting resistant to bursting susceptible cells for each experiment (means of corresponding quadruplicates). All LB agar plates and soft agar used throughout this study was supplemented with 25 μg/ml chloramphenicol. These CRISPR-Cas system efficiencies at different MOIs were tested if they are statistically different from each other using two-tailed unequal variances t-test at 0.05 confidence level using RStudio (R Core Team, 2013).

Determining the mean number of phages per cell

Request a detailed protocol

The cultures that were plated using standard plaque assays in the “Efficiency of the CRISPR-Cas system” experiment were also plated on LB agar plates containing 25 μg/ml chloramphenicol to determine the number of surviving CFUs. The numbers of bursting and surviving susceptible cells were subsequently used to determine the actual mean number of adsorbed phages per cell. The fraction of susceptible cells surviving the phage challenge experiment was assumed to correspond to the Poisson probability that a cell does not encounter any phage, which was than used to determine the mean of the Poisson distribution, which corresponds to the mean number of phages per cell.

Herd immunity in a liquid culture

Request a detailed protocol

Herd immunity in a liquid culture was tested in 200 μl of LB broth supplemented with 25 μg/ml chloramphenicol in Nunclon flat bottom 96 well plate in a Bio-Tek Synergy H1 Plate reader. Bacterial cultures were diluted 1 in 1000 and mixed in the following ratios of resistant to susceptible cells: 50:50, 75:25, 87.5:12.5, 93.75:6.25, 96.88:3.13, 98.44:1.56, 99.22:0.78, 99.61:0.39, 99.8:0.2, 99.9:0.1, 99.95:0.05, 100:0 %. T7 phage was added at a multiplicity of infection (MOI) of 10-4 (50 plaque forming units (pfu) per culture) to resemble an epidemic initiated by the burst size from one infected cell and the cultures were monitored at an optical density 600 nm for 18 hours post inoculation (hpi). Each population composition was replicated 16 times. Herd immunity in diluted LB was measured in LB broth mixed with 1X M9 salts in ratios 1:1 (50% LB) and 1:4 (20% LB) using the same protocol as for 100% LB broth. Each population composition was replicated 18 times.

Time-lapse imaging of plaque growth

Request a detailed protocol

Soft LB agar (0.7%) containing 25 μg/ml chloramphenicol was melted and 3 ml were poured into glass test tubes heated to 43C in a heating block. After the temperature equilibrated, 0.9 ml of a bacterial culture consisting of resistant and susceptible cells (ratios 10% – 100% of susceptible cells, 10% increments) were diluted 1 in 10 and added to the tubes. Then, 100 μl of bacteriophage stock, diluted such that there would be 10 plaques per plate, was added to the solution. Tubes were vortexed thoroughly and poured as an overlay on LB agar plates containing 25 μg/ml chloramphenicol. The plates were placed on scanners (Epson Perfection V600 Photo Scanner) and scanned every 20 minutes in ”Wide Transparency mode“ for 48 hours in 30C. A total of 3 scanners were employed with a total of 12 plates, plus a no phage control plate and 100% resistant control outside the scanners (see Appendix figure 3). No plaques were detected in the 100% resistant controls. Time-lapse images were used to calculate the increase of individual plaque areas using image analysis software PerkinElmer Volocity v6.3 and Fiji v1.0 (Schindelin et al., 2012).

Bacterial growth on soft agar

Request a detailed protocol

Growth rate of susceptible bacteria in soft LB agar (0.7%) was measured by sampling from a petri dish with a soft agar overlay with bacteria prepared in the same way as the plaque assays except an absence of the phage. Sampling was performed in spatially randomized quadruplicates at the beginning of the experiment and subsequently after 2, 4, 6, 8, 10, 12, 14, 16, 24, 32, 40, and 48 hours using sterile glass Pasteur pipettes (Fisherbrand art.no.: FB50251). Samples were blown out from the Pasteur pipette using an Accu-jet pro pipettor into 1 ml of M9 buffer pre-warmed to 43C, vortexed for 15 seconds and incubated for 10 minutes in 43C with two more vortexing steps after 5 and 10 minutes of incubation. Samples were serially diluted and plated on LB agar plates containing 25 μg/ml chloramphenicol. How bacterial densities change over time, measured as CFU/ml, is shown in Figure 4A.

Bacterial growth rates in liquid culture

Request a detailed protocol

Nutrient-dependent growth rate of susceptible bacteria was measured in Nunclon flat bottom 96 well plate in Bio-Tek Synergy H1 Plate reader in 30C. Overnight LB cultures were diluted 1:200 in media consisting of LB broth mixed with 1X M9 salts in ratios 10:90, 20:80, 30:70, 40:60, 50:50, 60:40, 70:30, 80:20, 90:10 and 100:0. Final volume was 200 μl. Optical density at 600 nm was measured every 10 min. Every treatment was replicated eight times. Natural logarithm of the optical density values was calculated to determine the growth rate using a maximal slope of a linear regression of a sliding window spanning 90 min.

The resulting growth rates for various nutrient concentrations fit well with Monod’s growth kinetics,

(9) α=αmaxNKc+N.

Results for the two fitting parameters, αmax and Kc, are listed in Table 1. The whole dataset, including the fit, is displayed in Figure 4B.

Test for a difference in growth rates of resistant and susceptible bacteria was done in LB broth in the same manner as nutrient-dependent growth rate measurements. Two-sample t-test was performed on acquired growth rate data at 0.05 confidence level using RStudio (R Core Team, 2013).

All growth media used in growth rate measurements were supplemented with 25 μg/ml chloramphenicol.

Phage burst sizes

Request a detailed protocol

Phage burst sizes in bacteria growing at different growth rates were measured by one-step phage growth experiments in LB mixed with 1X M9 salts in the following ratios 0:100, 20:80, 50:50 and 100:0. The burst sizes were calculated as the ratio of average number of plaques before burst to average number of plaques after burst. Consecutive samplings before and after burst were used for the calculation if they were not significantly different from each other (two sided t-test, p>0.05). All experiments were performed in triplicates.

Phage latent periods

Request a detailed protocol

Phage latent periods were determined from the phage burst size experiments as the time interval between the first and the last significantly different consecutive sampling between those used for phage burst size calculations.

Speed of phage expansion

Request a detailed protocol

The speed of the phage expansion was measured as difference in radii of plaques over time. Statistical tests allowed to infer that the reduction of expansion speed is significant already for small deviations from the 100% susceptible control experiment, as described and shown in Figure 7.

Phage diffusion in soft agar

Request a detailed protocol

Soft M9 salts soft agar (0.5%) was supplemented with SYBR safe staining (final conc. 1%) and poured into glass cuvettes (VWR type 6040-OG) to fill 2cm of the cuvette height. After soft agar solidification, the same stained soft agar was supplemented with T7 phage particles to a final concentration 1011pfu/ml and poured on top of the agar without phages. The cuvettes were monitored in 30C every hour for 40 hours at the SYBR safe emission spectrum peak wave length 524nm illuminated with the SYBR safe excitation spectrum peak wave length 509nm. The diffusion constant was estimated as the best fit parameter for the spread of fluorescent phages through the soft agar over time.

First we computed the luminosity Li of fluorescence (a gray-scale value defined as L=0.2126R+0.7152G+0.0722B from the RGB image) as average over the width of the cuvette for pixel row i, and corrected the profiles of luminosity Li by subtracting the background value. This background value was estimated as a linear fit at the end of the profile without phages, where only the gray value of the agar was measured. Moreover, in our setup luminosity seems to saturate at values above 0.4 where it does not have a simple linear dependence on fluorescence anymore: diffusion would lead to a decrease of the signal behind the inflection point of the profile and increase after the inflection point, but images only show increasing profiles – the bulk density does not decay. Thus, any estimate should only take the part of the profile that is below the threshold value of 0.4 into account (see Figure 8).

The diffusion constant D itself was estimated as the minimal value of the total squared deviation of the convoluted profile L(t) (at time t) with a heat kernel K(D) compared to the profile L(t+1) at time t+1,

(10) D=minDi((je-(i-j)2/4D4πDLj(t))-Li(t+1))2.

Such a convolution with the heat kernel Kij(D)=(4πD)-1/2exp(-(i-j)2/4D) assumes that the only change in the profile is due to diffusion for a time span of length 1 with i and j indices of pixels. Thus, expression Equation 10 estimates the diffusion constant in units of pixel2/frame, where frame is the time difference between two images. Several estimates are averaged over different snapshots in the whole experiment that spans 40h in intervals of 1h each.

The final estimate in appropriate units is

(11) D1.17(±0.26)102 mm2/h,

which is in agreement with previous measures of phage diffusion (Stent and Wollman, 1952; Bayer and DeBlois, 1974; Briandet et al., 2008).


Phage growth

Request a detailed protocol

In the main text we stated that relevant processes for phages growing on bacteria are given by the set of reactions Equation 2. In the following, we will analyze an extended version of our model, which takes all these processes into account. We try to justify our approximations and explain the reasoning behind leaving parts of the full model out. While reactions for single bacteria or phages are inherently stochastic in nature, we assume that the involved numbers are large enough such that the dynamics can be described with deterministic differential equations for the populations. Furthermore, reaction rates are identified with the inverse of the average time scale of the process. Thus, the full model is given by the coupled differential equations,

(12a)tBs =αBsA[Bs,P|Bs,Br],(12b)tBr =αBrA[Br,P|Bs,Br]+ρIr,(12c)tIs =A[Bs,P|Bs,Br](1/λ)Is,(12d)tIr =A[Br,P|Bs,Br](1/λ)IrρIr,(12e)tP =(β/λ)(Is+Ir)i{s,r}A[Bi,P|Bs,Br]i{s,r}A[Ii,P|Is,Ir],(12f)tN =α/y(Bs+Br).

Both bacterial populations Bi, i{s,r}, grow with rate α and decay via adsorption of phages A[Bi,P|Bs,Br], an expression that is specified below. Infected populations Ii gain numbers by adsorption and decrease via bursting. Resistant bacteria also can recover from their infected state with a recovery rate ρ. Phages grow by bursting cells, and lose numbers by adsorption to the various bacterial populations. Moreover, explicit dynamics for nutrients is considered, which are drained by each grown cell inversely proportional to the yield Y, the conversion factor between nutrient concentration and cell numbers. Essentially, this last equation acts as a timer, when we switch from abundant resources to the depleted state: all growth parameters change significantly upon nutrient depletion. Nevertheless, despite the possible deviations, we assume depletion time is given by the simple estimate Equation 3 and only treat the two possible states of abundant and depleted nutrients.

Adsorption of phages, given by the term A[Bi,P|Bs,Br], can be influenced by the whole distribution of populations within the culture. In liquid medium, a common assumption is that this term is proportional to the concentrations of both the phages and cells (Weitz, 2016),

(13) A[Bi,P|Bs,Br]=δBiP,

with an adsorption constant δ. This expression assumes constant mixing of the population and relatively short contact times between phages and bacteria. In general, this system of equations is akin to Lotka-Volterra dynamics, which has been analyzed in great detail, eg. (Hofbauer and Sigmund, 1998; Nowak, 2006).

For our ensuing analysis, we neglect the population of infected resistant bacteria Ir. Upon examining Equation 12 we find that most cells to leave their infected state by reducing phage DNA via CRISPR/Cas instead of bursting if ρ1/λ. If furthermore ρδP, which is true at least in the initial stages of the experiment, essentially all infected resistant bacteria immediately recover from a phage infection. Consequently, with both conditions, the resistant infected bacteria tend to vanish, Ir0, and their dynamics can be neglected. Only in the Appendix (section Simulation of recovery rate) we release this assumption to explicitly cover the full dynamics of Equation 12 in simulations to estimate values for ρ.

Exponentially growing bacteria lead to double exponential phage growth

Request a detailed protocol

For convenience, we transform the populations to the total bacterial density B=Bs+Br and introduce the fraction of susceptible cells S=Bs/B. The crucial assumption for the remainder of this section is that phages burst immediately after infection, λ=0, such that we can ignore all infected populations. While not a very biological condition, it allows to analyze the model in more detail. Using these simplifications, we obtain

(14a)tB =(αδSP)B,(14b)tS =S(1S)δP,(14c)tP =(βS1)δBP.

If we assume that in initial stages of phage growth the number of phages is small, ie. δPα𝒪(1h-1), the dynamics of bacteria and the fraction of susceptibles simplify to tB=αB and tS=0. Note that this term δP also occurs in the linear phage dynamics, where it cannot be neglected. In this instance, we need to view δB as a coefficient, which is likely much larger initially. This set of simplified equations can be solved in closed form,

(15a)S(t) =S0,(15b)B(t) =B0exp(αt),(15c)P(t) =P0exp((S0β1)δB0(exp(αt)1)/α).

The structure of phage dynamics is particularly important here – it exhibits a double-exponential dependence on time t, which is a very fast, almost explosive, growth. Such double-exponential growth leads to very large population sizes within a short amount of time (but after an extended initial delay). This general behavior of the solution is independent of the actual growth rate of phages, which only has to be positive. Thus, inspecting the exponent in Equation 15 yields the condition

(16) βS0>1

for phage growth to be positive. Incidentally, relation Equation 16 is the naive estimate for the number of successful additional infections arising from a single burst. The double exponential time-dependence is central for our arguing that the dynamics can be described by threshold phenomena, given by conditions like Equation 16: Usually, phages are negligible in the dynamics until they grow fast enough to large enough size, such that it is too late for the bacterial population to deal with the overwhelming phage population.

An important question in the context of these solutions is whether nutrients run out before this double-exponential growth of phages occurs. Hence, we compute the time Tδ defined as when phages reach a population of P(Tδ)=1/δ assuming phages grow as Equation 15 until then. After Tδ the assumptions that allowed to obtain Equation 15 are not valid anymore. Inverting Equation 15 for time leads to

(17) Tδ=1αlog(1+αlog(1/δP0)(βS0-1)δB0).

Subsequently, we can compare this estimate Tδ to the depletion time Tdepl=(1/α)log(B/B0). When rearranging the inequality Tdepl>Tδ in terms of the (initial) fraction of susceptibles S0, we obtain

(18) βS0>1+αlog(1/δP0)δ(BB0).

This expression Equation 18 is a condition for phages to reach “large” population sizes before nutrients are depleted by bacteria. The final population density B usually fulfills δB/α1, such that the correction given by the second term of Equation 18 can be considered small. Thus, if phages grow (βS0>1), they also grow very fast with a double-exponential time-dependence and reach a considerably large population size before bacteria stop multiplying (for almost all parameter values).

Extending analysis to finite burst times

Request a detailed protocol

The analysis above only treated the case λ0. However, we reported that the latency time λ increases significantly when bacterial growth rate α declines, see Table 2. Considering finite latency times entails dealing with an infected bacterial population I. (However, we identify IIs and set Ir=0.)

To this end, note that we can rearrange Equation 12 to (1+λt)I=λδSBP using the adsorption model in Equation 13. Hence, we can use the differential operator (1+λt) and apply it directly to Equation 12e to reduce the dependence on I in this equation at the cost of introducing higher order derivatives. In particular, we obtain

(19) λt2P+(1+λδB)tP+δB(βS-1-λα)P=0,

where we also inserted tBαB in the last term, as we aim again for a solution at initial times where δPα. The effects of the limit λ0 are directly observable – no terms are undefined in this limit. In particular, we find that equation Equation 19 and λ=0 lead directly to the dynamics of phages we just analyzed above, obtaining solution Equation 15.

In principle, Equation 19 is a hyperbolic reaction-diffusion-equation, which is known to occur upon transformation (or approximation) of time-delayed differential equations (Fort and Méndez, 2002b). For initial times we can use the solutions B(t)=B0exp(αt) and S(t)=S0. To proceed, we introduce the auxiliary variable

(20) z(t)=-δB0exp(αt)/α,

and assume P(z) as a function of this new variable z. We need to transform the differential operators of time derivatives, and obtain t=z(t)tz=αzz and t2=(αzz)(αzz)=α2(zz+z2z2). Inserting these expressions in Equation 19 and multiplying the whole equation with (α2λz)-1 yields the dynamics for phages,

(21) 0=zz2P(z)+(b-z)zP(z)-aP(z),

where the two extant constants are a=1-(βS0-1)/(λα) and b=1+1/(λα). Equation 21 is called 'Kummer’s equation' with confluent hypergeometric functions F11 as solutions (Abramowitz and Stegun, 1964, pg. 504),

(22) P(z)=AF11(a,b;z)+Bz1-bF11(a-b+1,2-b;z).

The two integration constants A and B can be determined via the initial conditions P(t=0)=P0 and (tP)(t=0)=-δB0P0. Using these conditions, the shape of the solution is again similar to before with λ=0 (double exponential time-dependence), although λ>0 introduces some skew. The most important aspect of this solution Equation 22 is to compute the parameter combination where it switches from a decreasing to increasing function over time. A careful analysis reveals that at the parameter value a=0 the behavior of the solution changes. Consequently, we find the condition for growing phage populations,

(23) βS0>1+λα,

which is a non-trivial extension including finite latency times λ.

Note, however, that this relation Equation 23 does not indicate a correction to the general epidemiological parameter R0, which can be identified with β in our model. Rather, it shows that a growing bacterial population requires the phage population to grow even faster for a continuous chain of infections in an epidemic. The term λα denotes the ratio of generation times of pathogen over host, which in most cases is small and negligible compared to 1. For bacteria and phages, however, which have similar generation times, such a correction is needed to describe the effects of growing host population sizes. In contrast, many other epidemiological models assume the host population size constant and only pathogens are increasing (or decreasing) in number.

While our result Equation 23 suggest that it also should hold in the limit α0, it might not necessarily be so. This specific limit is actually quite important for the time when nutrients are depleted in the experiments. However, at several instances in the calculations above we implied a positive α>0. The most important of these is the transformation to z(t)=-δB(t)/α, which actually exhibits two problems: dividing by α should not be allowed and B(t) is essentially constant and cannot serve as a variable in a differential equation. We also neglected the second term in tB=(α-δSP)B throughout our calculation. For α=0 this second term is dominant in bacterial dynamics and would generate non-linear phage dynamics if inserted for tB right before stating Equation 19. However, we expect that albeit the process will run very slow, and might not be measurable in experiments, the simple condition βS0>1 could indicate phage expansion and bacterial decay.

Growth of phages on plated bacterial lawn

Request a detailed protocol

Spatial modelling of phage expansion has produced several predictions for how plaque radius r and expansion speed v are influenced by experimentally adjustable parameters (Kaplan et al., 1981; Yin and McCaskill, 1992; You and Yin, 1999; Fort and Méndez, 2002a; Ortega-Cejas et al., 2004; Abedon and Culler, 2007; Mitarai et al., 2016). Here, we try to use a minimal model to estimate these two observables, based on the considerations of previous sections.

One of the main complications arises from the fact that all densities in Equation 12 have a spatial dimension in addition to their time dependence, Bi=Bi(x,t), i{s,r}. As explained in the main text we only consider phage diffusion, heterogeneities in all other densities are generated only by phage growth. The additional spatial dimension imposes a particular contact network between phages and bacteria, which are not entirely random encounters anymore: One can expect that the size of the bacterial neighborhood B^ phages are able to explore is only slightly determined by the actual density B, and can be assumed constant for most of the experiment, B^(B)const. Consequently, the adsorption term can be written in the following way,

(24) A[Bi,P|Bs,Br]=δBiBs+BrP,i{s,r},

which only depends on the relative frequencies of bacterial strains. The adsorption constant δ is both the rate of adsorption and inter-host transit time as determined by the diffusion constant D. Thus, one can expect the formal dependence δ=δ(D,B^(B)). For our particular experimental setup, however, δ will be treated as a constant. This adsorption term Equation 24 leads to the dynamics of phages

(25) tP=D2P+G[P,S],

where we collected all contributions to phage growth in G[P,S] and added the spatial diffusion term D2P. For simplicity, we consider only expansion in a single dimension (2x2), which has been found to coincide well with the dynamics of plaque growth (Yin and McCaskill, 1992). The growth term for phages is then defined as,

(26) G[P,S]=δ(Sβ-1-λα)P,

where we also consider the correction λα obtained from the analysis in liquid culture. Due to the different absorption dynamics on plates, however, this correction might be a slight overestimate of the actual term that accounts for bacterial growth. Reaction-diffusion equations similar to Equation 25 have been first analyzed about 80 years ago (Fisher, 1937; Kolmogorov et al., 1991) and since then treated extensively, e.g. (Murray, 2002; van Saarloos, 2003). They admit a traveling wave solution – here, this corresponds to phages sweeping over an uninfected bacterial lawn. In general, the asymptotic expansion speed for the traveling wave solutions is given by the following expression,

(27) v=2D(PG)[0,S]=2DδSβ1λα.

Only the linearized growth rate of phages at very low densities is relevant for the expansion speed, PG[P=0,S]. Thus, the fraction of susceptible individuals S should be unchanged from its initial value S0. It should be noted, that only for S0β>1+λα does Equation 27 remain valid, otherwise we have v=0. Such a scenario is relevant when nutrients are depleted and phage growth parameters changes to βdepl and λdepl.

The expression for the expansion speed also shows the need for the spatial adsorption model in Equation 24, in contrast to the liquid case Equation 13. If adsorption would directly depend on the bacterial density B, the additional linear dependence on B in Equation 26 would lead to an exponentially increasing speed during the experiment. This is in clear contradiction to experimental observations.

The density of phages behind the expanding front is large and as previously noted at large MOIs the CRISPR-Cas system fails to provide effective immunity (see section Materials and methods and appendix Infection load and efficiency of the CRISPR/Cas system). However, in comparison to an un-structured environment (e.g., liquid) the structured environment effectively limits transit of phage from within a plaque to the expanding front: The combined effect of growth and diffusion usually generates a much faster expansion of phages during plaque formation, than diffusion alone. Only when nutrients are depleted, can pure diffusion processes explain the slow decrease in speed observed in experiments (see Figure 7A). Our model assumes a sharp drop to v=0 at Tdepl for small S.

In order to derive an expression for the plaque radius r, we integrate the expansion speed Equation 7 over time, r(t)=0tdtv(t). Employing the simplification that only two values of phage growth are necessary to describe the dynamics – before Tdepl phages grow normally with β and λ, after Tdepl phage growth changes to βdepl and λdepl – we can evaluate the integral for the radius directly, arriving at,

(28) r(t)={2tDδSβ1λα,0<t<Tdepl,2Dδ(TdeplSβ1λα+(tTdepl)Sβdepl1),Tdepl<t.

Using this expression we estimated the adsorption constant δ from the growth experiments as it difficult to measure in practice. This estimate is done for radii exactly at the time of nutrient depletion Tdepl, and excluding the control experiment with only susceptible cells.

Predictions of our model show a discrepancy from experimental results on plates after depletion. We independently estimated βdepl=3.0, which results in Hdepl=(βdepl-1)/βdepl0.67. Thus, all experiments with S>0.33 should exhibit expanding plaques after nutrients are depleted. In the experimental setup plaques stop expanding in all mixtures of resistant to susceptible cells (S0.9), which would correspond to βdepl<1.1. This value is, however, still within experimental accuracy of our estimates of βdepl.

Appendix 1

Additional theoretical considerations

Simulation of recovery rate

Throughout the main text we assumed that resistant bacteria are completely immune to phage infection as their CRISPR/Cas system immediately kills adsorbed phages. However, experimental observation suggest that for fractions close to what we predicted as herd immunity threshold, all bacteria eventually die. Thus, in the following section we use numerical simulations to investigate the full set of equations Equation 12, with a particular focus on the question why the whole bacterial population goes extinct. As it turns out, this requires using finite values for the recovery rate ρ (instead of the ρ approximation employed previously).

A major difficulty in analyzing the full model Equation 12 is finding appropriate parameter values. In particular, we need values for the adsorption constant δ, the recovery rate ρ and the yield coefficient Y. Undiluted LB medium is known to support a population of 5109cells/ml. Thus one can easily estimate Y as the inverse of this number, when nutrients are measured in units of dilutions, which we already used throughout this publication (undiluted medium corresponds to N=1). Parameter scans in simulations reveal that the actual value of the adsorption constant δ does usually not influence the actual outcome (collapsed or surviving bacterial population), it only adjusts time scales. However, deviations in time scales are insignificant, even when δ is changed by orders of magnitude, δ𝒪(10-610-8). They are roughly an hour or less, which is small compared to the expected duration of the experiment that lasts a few hours. For definiteness, we use the value of δ=10-7h-1 for our simulations. That the value of the adsorption constant has only a minor impact on phage growth on bacterial cultures, is also in line with previous findings (Mitarai et al., 2016).

The most elusive parameter is the recovery rate ρ. A first indication of the value of ρ can be drawn from our experiments on bursting resistant cells, summarized in Figure 2. As the probability for bursting resistant cells is 3 orders of magnitude smaller than for susceptible bacteria, we can use 1/λ𝒪(1) to estimate ρ𝒪(103). However, our results also indicate that recovery via the CRISPR/Cas system heavily depends on MOI, implying that ρ depends on the actual densities of phages and bacteria. Nevertheless, as experimental determination of recovery is complicated, even more so determining a functional dependence on dynamically changing densities B and P, we assume that ρ is constant.

We ran parameter sweeps in simulations and compared the outcome – collapsed or surviving bacterial populations – to the observed experimental results (see Figure 3). The best agreement of simulations and experiments was reached with ρ𝒪(1). Lower values of ρ do not allow the resistant population to recover from phage infection, while for larger values of ρ, phages are drained from the culture very fast. Such a small value of ρ is most likely related to the recovery at very large MOI, when the densities involved in the dynamics are large, which dominate the overall observed dynamics. At this time phages repeatedly infect the same bacteria and their CRISPR/Cas immune system cannot deal with such an infection load (or only too slow). Thus, we can argue that our final choice ρ=1.5h-1 is the recovery rate when the CRISPR/Cas system is heavily stressed, which is comparable to the actual burst rate 1/λ for phages.

In Appendix 1—figure 1 we show three exemplary sets of trajectories for bacteria and phage. For a tiny fraction of susceptibles, S=10-3, which is well below the herd immunity threshold (see Figure 3), phages do not thrive on the limited number of favorable hosts and decay fast after a slight increase initially. For intermediate fractions of susceptibles, S=0.04, we observe more complex, non-monotonic trajectories of bacterial populations. For such values of S we also observe mixed outcomes in experiments, see Figure 3. When S is increased further (S=0.06), enough susceptible bacteria exists to produce enough phages and ultimately the whole bacterial population goes extinct.

Appendix 1—figure 1
Simulated trajectories for all populations in liquid culture for the extended model, including infected and recovering bacteria.

Trajectories are obtained by numerically integrating equations Equation 12, using parameters listed in Appendix 1—table 1 and additionally N=1, Y=210-10cells-1, δ=10-7h-1 and ρ=1.5h-1. (A) For population compositions with a large majority of resistant cells (S=10-3), phages get wiped out fast. (B) For intermediate S (close to parameters where we observe both, collapsed and surviving, populations, see Figure 3), the populations exhibit a complex, non-monotonic trajectory. After fast initial growth of phages, bacterial populations decay but ultimately can recover. (C) If the fraction of susceptibles is too large (S=0.06), the whole bacterial population is infected and succumbs to the overwhelming phage infection. See supporting text for more detailed information.

Appendix 1—table 1
Parameters used in simulations shown in Appendix 1—figure 1.
Bacterial growth ratea0.631/hTable 1, Figure 4
YieldY210-11/cellMeasured in dilution of LB
Recovery rateρ1.51/hSee this appendix
Adsorption constantδ10-71/(hcell)See this appendix
Diffusion constantD1.1710-2mm2/hSee Methods
Burst sizeβ85.6phages/cellTable 2, Figure 5
Latency timeλ0.60hTable 2, Figure 5
Initial bacterial populationB0105cells
Initial phage populationP010phages

The purpose of the extended model in this section was to justify the fact that phages can wipe out the whole bacterial population, which was not possible in the simplified model used in the main text. There, the resistant bacterial population was basically unaffected by phages and just acted as ‘sink’ for phages. However, also in this extended model, we see a very similar behavior in terms of the threshold phenomena reported earlier in the manuscript.

Infection load and efficiency of the CRISPR/Cas system

In the section Modelling we showed that positive phage growth leads eventually to a very fast increase in the phage population, that occurs before nutrients are depleted (for almost all realistic parameters). This behavior of the dynamics was also observed in the extended simulation model presented in the last section. Moreover, as a condition we used that the phage population reaches a size P1/δ, which is after all arbitrary – it only determines if we can employ useful simplifications and approximations to model equations. However, simulation results presented in the last section indicate that the bacterial population starts to decay soon after such a threshold P1/δ is exceeded.

In order to proceed, we investigate the system at time Tδ further. We assume that the phage population is large enough that it will not be degraded by the CRISPR/Cas immune system. The threat to immediate phage extinction is low at this point. The actual equations are hard to solve directly, hence we revert to simple balance equations, ignoring the dynamical component. Specifically, we compare the number of (present and eventually produced) phages to the number of infections needed to wipe out the whole population. To incorporate the effects of the bacterial immune system in resistant bacteria, we assume that they need M>1 infections before they burst and produce only κβ phages, which reduces the burst size by a (yet unspecified) factor 0<κ<1. κ=1 implies that resistant cells produce the same number of phages as susceptible cells, while κ=0 indicates only cell death. Combining these considerations yields

(29) 1/δphages present+βS0B(Tδ)phage production Bs+κβ(1S0)B(Tδ)phage production Br>S0B(Tδ)infections Bs+M(1S0)B(Tδ)infections Br,

where the left side indicates the total number of phages, while the right side indicates the number of necessary infections to kill all bacteria. The number of bacteria B(Tδ) can be estimated by inserting the time Tδ from Equation 18 into the exponential growth Equation 15b. Subsequently, we can rearrange Equation 29, obtaining a bound on M:

(30) M<1/δB(Tδ)+S0(β1)1S0+κβ.

The first term 1/δB(Tδ) indicates the ratio of phages to bacteria at time Tδ, and can be considered small for non-extremal parameters compared to the other terms. This fact justifies our assumption that the actual value of δ is not crucial. This number M might allow some insight into the effectiveness of the CRISPR/Cas immune system. For a fraction of susceptibles S=0.03, which corresponds to the minimal value where we observe only collapsed bacterial populations in undiluted LB medium (see Figures 3,5), we would obtain the relation M3+86κ. Thus, each resistant bacterial cell could degrade up to 𝒪(101102) phages before their CRISPR/Cas system cannot cope with the infection load anymore.

Appendix 2

Supplementary results

Reduction in number of plaques in spatially structured populations

The reduction in the number of plaques with increasing proportions of resistant bacteria is shown in Appendix 2—figure 1.

Appendix 2—figure 1
Number of plaques declines faster than proportionally to the fraction of resistant bacteria.

Number of plaque forming units observed on a plate (y-axis) for different proportions of resistant bacteria (x-axis). Grey numbers below each boxplot indicate the average number of phages inoculated in the respective treatment. The numbers of phages inoculated were chosen to retain the expected number of plaques on the plate (green dashed line) as in the 0% resistant treatment (red boxplot). Plates were prepared using identical procedure as in Time-lapse imaging of plaque growth (see Materials and methods). The data presented in this figure can be found in Appendix 2—figure 1—source data 1.

Appendix 2—figure 1—source data 1

Measurements of plaque numbers in populations consisting of varying proportions of resistant to susceptible bacteria.


Appendix 3

Additional information on experimental setup

Our experimental setup for the scanner system is shown in Appendix 3—figure 1.

Appendix 3—figure 1
Image of the scanner system.

Photograph of the scanner system used for time-lapse imaging of phage spread in spatially structured bacterial populations. Three scanners (Epson Perfection V600 Photo Scanner) simultaneously scanned 12 plates in total every 20 min in 30C for 48 hr per experiment.


Appendix 4

Source data and code

Time-lapse images of spread of T7 phage epidemics in Escherichia coli spatially structured populations are available on the Dryad Digital Repository: https://doi.org/10.5061/dryad.42n44. Source code of the model presented here is available on GitHub (https://github.com/lukasgeyrhofer/phagegrowth) (Payne et al., 2018) and its archived version is accessible through Zenodo: https://dx.doi.org/10.5281/zenodo.1038582 and https://github.com/elifesciences-publications/phagegrowth.

Data availability

The following data sets were generated


  1. Book
    1. Abramowitz M
    2. Stegun IA
    Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, 55
    Courier Corporation.
  2. Book
    1. Anderson RM
    2. May RM
    (August 1992)
    Infectious Diseases of Humans: Dynamics and Control
    OUP Oxford.
    1. Bayer ME
    2. DeBlois RW
    Diffusion constant and dimension of bacteriophage phi X174 as determined by self-beat laser light spectroscopy and electron microscopy
    Journal of Virology 14:975–980.
    1. Blumenthal RM
    2. Cheng X
    Modern Bicrobial Genetics
    177–225, Restriction-modification systems, Modern Bicrobial Genetics, 10.1002/047122197X.ch7.
    1. Fenner F
    Smallpox: emergence, global spread, and eradication
    History and Philosophy of the Life Sciences 15:397–420.
  3. Book
    1. Hamer WH
    Epidemic Disease in England: The Evidence of Variability and of Persistency of Type
    Bedford Press.
  4. Book
    1. Kolmogorov A
    2. Petrovsky I
    3. Piscounoff N
    A Study of the Diffusion Equation with Growth of the Quantity of Matter and Its Application to a Biology Problem
    In: Tikhomirov VM, editors. Selected Works of A. N. Kolmogorov, 1. Dordrecht, The Netherlands: Kluwer Academic Publishers. pp. 242–270.
    1. Murray JD
    2. Stanley EA
    3. Brown DL
    (1986) On the spatial spread of rabies among foxes
    Proceedings of the Royal Society B: Biological Sciences 229:111–150.
  5. Book
    1. Murray JD
    Mathematical Biology I: An Introduction, Vol. 17 of Interdisciplinary Applied Mathematics
    New York: Springer.
    1. Nordström K
    2. Forsgren A
    Effect of protein A on adsorption of bacteriophages to Staphylococcus aureus
    Journal of Virology 14:198–202.
  6. Book
    1. Nowak MA
    Evolutionary Dynamics
    Harvard University Press.
    1. Payne P
    2. Geyrhofer L
    3. Barton NH
    4. Bollback JP
    (2018) phagegrowth
    phagegrowth, 42ec97e, https://github.com/lukasgeyrhofer/phagegrowth.
  7. Software
    1. R Core Team
    (2013) R: a language and environment for statistical computing
    R Foundation for Statistical Computing, Vienna, Austria.
  8. Book
    1. Weitz JS
    Quantitative Viral Ecology: Dynamics of Viruses and Their Microbial Hosts
    Princeton, United States: Princeton University Press.

Decision letter

  1. Richard A Neher
    Reviewing Editor; University of Basel, Switzerland

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "CRISPR-based herd immunity limits phage epidemics in bacterial populations" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom Richard A Neher is a member of our Board of Reviewing Editors and the evaluation has been overseen by and Arup Chakraborty as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Timothy Cooper (Reviewer #2).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.


This manuscript presents a mathematical model and a set of experiments examining how CRISPR-Cas mediated immunity to phage infection can result in herd immunity preventing the spread of the pathogen. The experiments are performed both in well-mixed liquid culture as well as in structured environments on agar. The critical herd immunity thresholds and the plaque radii follow model predictions as the parameters of the system as varied. All reviewers appreciated the elegant combination of experiments and theory, but we would like the authors to address the following points:

Essential revisions:

1) The system is set up such that resistant cells don't grow once infected. We would like to see the importance of this assumption discussed. How would model predictions change if cells kept growing? How important is rapid degradation of phage by resistant cells?

2) We would like to encourage the authors to move some derivations presented in Materials and methods section to the Results section. While detailed algebra is better kept in the Materials and methods section, additional mathematical details would help readers appreciate the theoretical results without having to jump back and forth to the Materials and methods section.

3) The experiments are started with a small inoccula which result in a broad crossover at the herd immunity threshold. Is the observed stochasticity consistent with that expected from a small founder population? Or are there additional sources of stochasticity? Would the transition in Figure 4C be sharper if larger inoccula had been used? What is the nature of the error bars in Figure 3 and Figure 4C? Do they quantify uncertainty in the parameter measurements or the stochasticity of the dynamics?

4) There is a sophisticated body of theory on pathogen spreading with long range jumps (e.g. Hallatschek and Fisher, 2014). Can you provide some intuition/discussion how the spreading and herd immunity thresholds change as dispersal goes from 2D diffusion to occasional long-range jumps ultimately to the well-mixed case?

5) It is mentioned in the Abstract that herd immunity might facilitate coexistence between susceptible and resistant variants. This possibility is mentioned again in the discussion in one sentence but not elaborated on. While it seems obvious that selection for resistant variants is reduced when the pathogen can no longer spread due to herd immunity, it is not obvious that this herd immunity results in stable coexistence. Resistant variants would still sweep to high frequencies in areas with high pathogen load. How would a cost of resistance change the model behavior?

The quantitative model the authors parameterized with experiments could be used to investigate many of these points in greater detail. If the model for example predicted the width of the cross-over at the herd immunity threshold correctly, such comparisons would increase confidence in both the model and the experimental system. Furthermore, the model could be used to assess robustness to assumptions like growth arrest and rapid phage degradation.


Author response

Essential revisions:

1) The system is set up such that resistant cells don't grow once infected. We would like to see the importance of this assumption discussed. How would model predictions change if cells kept growing? How important is rapid degradation of phage by resistant cells?

We agree with the reviewers that this should have been discussed in more detail. Therefore, we added to the Results section a discussion of the importance of this property of our experimental system – namely, that resistant bacteria stop growing once infected.

Degradation of the phage by resistant bacteria is, in our opinion, an important pre-requisite of (bacterial) herd immunity. If the phage particles survive an encounter with a resistant bacterium (e.g. with a mutated receptor protein to which phages cannot adsorb) and can go on to infect other individuals, they would have to decay in the environment at the timescale comparable to their adsorption rate in order for resistant bacteria to provide any indirect protection to the susceptible ones. We argue that the usual persistence of the phage particles in both laboratory and natural environments must be substantially longer as dwelling outside the host likely covers a majority of phage life history. Also, common laboratory experience shows that titres of many different phage species remain stable over very large time periods in a variety of media and buffers. In order to clarify this point, we comment on the potential degradation of phage particles in the environment in the Introduction.

2) We would like to encourage the authors to move some derivations presented in Materials and methods section to the Results section. While detailed algebra is better kept in the Materials and methods section, additional mathematical details would help readers appreciate the theoretical results without having to jump back and forth to the Materials and methods section.

We have brought Equations 24 and 26 in a combined form, as well as their description, into the Results section to aid readers' understanding of the model. Both equations remain in the Materials and methods section to aid with understanding the derivations. In addition, we have added an explanation in the Results section of how the plaque radius was derived using the Equation 7 for the speed of expansion and Equation 3 for nutrient levels (subsection “Modelling herd immunity in spatially structured populations”). We also reformulated some text in the modelling section of the Results section and of the Materials and methods section for a better understandability.

3) The experiments are started with a small inoccula which result in a broad crossover at the herd immunity threshold. Is the observed stochasticity consistent with that expected from a small founder population? Or are there additional sources of stochasticity? Would the transition in Figure 4C be sharper if larger inoccula had been used? What is the nature of the error bars in Figure 3 and Figure 4C? Do they quantify uncertainty in the parameter measurements or the stochasticity of the dynamics?

We have now added a paragraph in the Results section explaining in detail the sources of stochasticity, the nature of the error bars in Figure 3 and Figure 4C and have commented on the potential effect of a larger phage inocula.

4) There is a sophisticated body of theory on pathogen spreading with long range jumps (e.g. Hallatschek and Fisher, 2014). Can you provide some intuition/discussion how the spreading and herd immunity thresholds change as dispersal goes from 2D diffusion to occasional long-range jumps ultimately to the well-mixed case?

In general, we would not expect that long range jumps change the herd immunity threshold H(α) itself. Spread of pathogens still stops when the fraction of susceptible hosts S is small such that βS < 1 + λα, and will continue as long as βS > 1 + λα is fulfilled. However, we expect that long-range jumps of phages will enhance the overall speed of expansion, and significantly increase the area covered by the phage infection. In order to address these concerns, we have added a paragraph (subsection “Modelling herd immunity in spatially structured populations”) providing some intuition of how long-range jumps might affect epidemic spread and the herd immunity threshold, in both general terms and in our experimental system.

5) It is mentioned in the Abstract that herd immunity might facilitate coexistence between susceptible and resistant variants. This possibility is mentioned again in the discussion in one sentence but not elaborated on. While it seems obvious that selection for resistant variants is reduced when the pathogen can no longer spread due to herd immunity, it is not obvious that this herd immunity results in stable coexistence. Resistant variants would still sweep to high frequencies in areas with high pathogen load. How would a cost of resistance change the model behavior?

We have elaborated on the role of herd immunity in coexistence of resistant and susceptible individuals as well as the pathogen in a population in greater detail in the Discussion section. We also discussed what role a cost of immunity may play in the dynamics of herd immunity and coexistence of the pathogen and the hosts.

The quantitative model the authors parameterized with experiments could be used to investigate many of these points in greater detail. If the model for example predicted the width of the cross-over at the herd immunity threshold correctly, such comparisons would increase confidence in both the model and the experimental system. Furthermore, the model could be used to assess robustness to assumptions like growth arrest and rapid phage degradation.


Article and author information

Author details

  1. Pavel Payne

    1. Institute of Integrative Biology, University of Liverpool, Liverpool, United Kingdom
    2. Institute of Science and Technology Austria, Klosterneuburg, Austria
    Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2711-9453
  2. Lukas Geyrhofer

    Department of Chemical Engineering, Technion - Israel Institute of Technology, Haifa, Israel
    Conceptualization, Data curation, Software, Formal analysis, Visualization, Methodology, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8043-2975
  3. Nicholas H Barton

    Institute of Science and Technology Austria, Klosterneuburg, Austria
    Conceptualization, Resources, Supervision, Methodology, Project administration, Writing—review and editing
    Competing interests
    No competing interests declared
  4. Jonathan P Bollback

    1. Institute of Integrative Biology, University of Liverpool, Liverpool, United Kingdom
    2. Institute of Science and Technology Austria, Klosterneuburg, Austria
    Conceptualization, Resources, Supervision, Funding acquisition, Validation, Methodology, Project administration, Writing—review and editing
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4624-4612


H2020 European Research Council (EVOLHGT No. 648440)

  • Jonathan P Bollback

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.


We are grateful to Remy Chait for his help and assistance with establishing our experimental setups and to Tobias Bergmiller for valuable insights into some specific experimental details. We thank Luciano Marraffini for donating us the pCas9 plasmid used in this study. We also want to express our gratitude to Seth Barribeau, Andrea Betancourt, Călin Guet, Mato Lagator, Tiago Paixão and Maroš Pleška for valuable discussions on the manuscript. Finally, we would like to thank the eLife editors and reviewers for their helpful comments and suggestions.

Reviewing Editor

  1. Richard A Neher, University of Basel, Switzerland

Publication history

  1. Received: October 7, 2017
  2. Accepted: March 8, 2018
  3. Accepted Manuscript published: March 9, 2018 (version 1)
  4. Version of Record published: April 27, 2018 (version 2)
  5. Version of Record updated: June 11, 2018 (version 3)


© 2018, Payne et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.


  • 2,827
    Page views
  • 424
  • 13

Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.

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)

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

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

Further reading

    1. Ecology
    Jakob Thyrring, Lloyd S Peck
    Research Article Updated

    Whether global latitudinal diversity gradients exist in rocky intertidal α-diversity and across functional groups remains unknown. Using literature data from 433 intertidal sites, we investigated α-diversity patterns across 155° of latitude, and whether local-scale or global-scale structuring processes control α-diversity. We, furthermore, investigated how the relative composition of functional groups changes with latitude. α-Diversity differed among hemispheres with a mid-latitudinal peak in the north, and a non-significant unimodal pattern in the south, but there was no support for a tropical-to-polar decrease in α-diversity. Although global-scale drivers had no discernible effect, the local-scale drivers significantly affected α-diversity, and our results reveal that latitudinal diversity gradients are outweighed by local processes. In contrast to α-diversity patterns, species richness of three functional groups (predators, grazers, and suspension feeders) declined with latitude, coinciding with an inverse gradient in algae. Polar and tropical intertidal data were sparse, and more sampling is required to improve knowledge of marine biodiversity.

    1. Ecology
    2. Evolutionary Biology
    Aspen T Reese et al.
    Research Article

    Domesticated animals experienced profound changes in diet, environment, and social interactions that likely shaped their gut microbiota and were potentially analogous to ecological changes experienced by humans during industrialization. Comparing the gut microbiota of wild and domesticated mammals plus chimpanzees and humans, we found a strong signal of domestication in overall gut microbial community composition and similar changes in composition with domestication and industrialization. Reciprocal diet switches within mouse and canid dyads demonstrated the critical role of diet in shaping the domesticated gut microbiota. Notably, we succeeded in recovering wild-like microbiota in domesticated mice through experimental colonization. Although fundamentally different processes, we conclude that domestication and industrialization have impacted the gut microbiota in related ways, likely through shared ecological change. Our findings highlight the utility, and limitations, of domesticated animal models for human research and the importance of studying wild animals and non-industrialized humans for interrogating signals of host–microbial coevolution.