Abstract
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.
https://doi.org/10.7554/eLife.32035.001eLife 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.
https://doi.org/10.7554/eLife.32035.002Introduction
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, ${R}_{0}$, 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 ${R}_{0}$ see [Heesterbeek, 2002]). Thus, ${R}_{0}$ indicates the epidemic spreading potential in a population. Given ${R}_{0}$ the herd immunity threshold is defined as,
which determines the required minimum fraction of resistant individuals needed to halt the spread of an epidemic. ${R}_{0}$ 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, 199208), 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 (DestoumieuxGarzón et al., 2005). Alongside these bacteria have evolved innate immune systems that target phage genomes for destruction. These include host restrictionmodification systems (RMS) (Blumenthal and Cheng, 2002), argonautebased RNAilike systems (Swarts et al., 2014), and bacteriophageexclusion (BREX) systems (Goldfarb et al., 2015). In addition to innate systems, bacteria have evolved an adaptive immune system called CRISPRCas (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, argonautebased RNAilike, and the CRISPRCas 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 CRISPRCas 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.
In addition to immune systemspecific 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 phagebacterial 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 CRISPRbased 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.
Results
Properties of resistant individuals
We engineered a resistant E. coli strain by introducing the CRISPRCas 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).

Figure 2—source data 1
 https://doi.org/10.7554/eLife.32035.005
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 timepoint 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).

Figure 3—source data 1
 https://doi.org/10.7554/eLife.32035.007
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%.

Figure 4—source data 1
 https://doi.org/10.7554/eLife.32035.009

Figure 5—source data 1
 https://doi.org/10.7554/eLife.32035.011

Figure 5—source data 2
 https://doi.org/10.7554/eLife.32035.012
From the experimental observations of the herd immunity threshold values we infer the phage ${R}_{0}$ using Equation 1. In an undiluted growth medium the phage ${R}_{0}$ 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 microcolonies, 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 realtime imaging (Figure 6A).

Figure 6—source data 1
 https://doi.org/10.7554/eLife.32035.015
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$, ${t}_{53}$ = $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$, ${t}_{32}$ = $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$, ${t}_{62}$ = $1.79$), and $5.67hpi$ with 30% ($p=0.0286$, ${t}_{53}$ = $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.

Figure 7—source data 1
 https://doi.org/10.7554/eLife.32035.017
The results presented in this and the previous section would allow us to use Equation 1 to infer a value for ${R}_{0}$ 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 ${R}_{0}$) 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 $\beta $ and latent period $\lambda $ depend strongly on the bacterial growth rate $\alpha $ (see Table 1).
The main processes in our model system can be defined by the following set of reactions,
Susceptible (${B}_{\mathrm{s}}$) and resistant (${B}_{\mathrm{r}}$) cells grow at a rate $\alpha $ (no significant difference in growth rate between strains, $\alpha ({B}_{\mathrm{s}})=0.551\pm 0.045{h}^{1}$, $\alpha ({B}_{\mathrm{r}})=0.535\pm 0.023{h}^{1},{t}_{70}=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 $\beta $ phage with a rate inversely proportional to the average latency $\lambda $. In contrast, resistant bacteria either survive by restricting phage DNA via their CRISPRCas 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\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{T}_{\mathrm{d}\mathrm{e}\mathrm{p}\mathrm{l}})=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,
Here, ${B}_{0}$ and ${B}_{\mathrm{\infty}}$ are the initial and final bacterial densities, respectively. In the initial exponential growth phase, our estimates from experimental data for growth parameters are $\alpha =0.63{h}^{1}$, $\beta =85.6\mathrm{phages}/\mathrm{cell}$ and $\lambda =0.60h$, for bacteria and phages, respectively (Tables 1 and 2). After time ${T}_{\mathrm{depl}}$, bacterial growth rate is set to zero ($\alpha =0$) and phage growth is reduced to ${\beta}_{\mathrm{depl}}=3.0\mathrm{phages}/\mathrm{cell}$ and ${\lambda}_{\mathrm{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 $\beta \phantom{\rule{thinmathspace}{0ex}}S\phantom{\rule{thinmathspace}{0ex}}>1$, which leads to a continuous chain of infections: the product of burst size $\beta $ 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 ${R}_{0}$ with the burst size $\beta $, 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 ${R}_{0}$ and $H$ as shown in Equation 1. We capture this dynamical effect in a correction to the previous estimate as $\beta S\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1+\lambda \alpha$ (see Materials and methods): more phages have to be produced for the chain of infections to persist in growing populations. The correction $\frac{\lambda}{1/\alpha}$ is the ratio of generation times of phages over bacteria – usually, such a correction is very small for nonmicrobial hosts and can be neglected. Ultimately, herd immunity is achieved if the threshold defined by $H=1{S}_{c}$ is exceeded, with ${S}_{c}$ the critical value in the inequality above. Rearranging, we obtain an expression for the herd immunity threshold
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 nonmonotonic trajectories for ${B}_{\mathrm{s}}$ and ${B}_{\mathrm{r}}$ (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 $\alpha $, phage burst size $\beta $ and phage latent period $\lambda $ 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, $\beta (\alpha )$, using a numerical quadratic fit (Figure 5A). Similarly, we estimated the dependence of the phage latent period on the bacterial growth rate, $\lambda (\alpha )$ (Figure 5B). Using these estimates we calculated the expected growth rate–dependent herd immunity threshold
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 $\beta $, $\lambda $ and $\alpha $. 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; OrtegaCejas 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\phantom{\rule{thinmathspace}{0ex}}(\pm \phantom{\rule{thinmathspace}{0ex}}0.26)\cdot {10}^{2}\text{}\phantom{\rule{thinmathspace}{0ex}}m{m}^{2}/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 ${\delta}^{\star}$.
Taking these considerations together, allows to write a reactiondiffusion dynamics for growth of phages $P$ on the growing bacterial population as
The first term accounts for the diffusive spread of phages, while the second term describes phage growth. This second term includes the correction $\lambda \alpha $ 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
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, timedependence of the speed only enters via the fraction of susceptibles $S$, which is assumed to stay at the initial ${S}_{0}$ 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 $r\approx v{T}_{\mathrm{depl}}$, with ${T}_{\mathrm{depl}}$ given by Equation 3. Using this expression we estimated the adsorption constant ${\delta}^{\star}$ 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 ${\delta}^{\star}=4.89(\pm 0.19)\cdot {10}^{2}\mathrm{bacteria}/\mathrm{phage}h$ 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,
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 pathogenhost system. Ultimately, this relative speed approaches zero with a universal exponent of $1/2$, when the fraction of resistant individuals $1S$ 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 ($\approx 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(\alpha )$ itself. Spread of pathogens still stops when the fraction of susceptible hosts $S$ is small such that $\beta S\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}1+\lambda \alpha$, and will continue as long as $\beta S\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1+\lambda \alpha$ is fulfilled.
Discussion
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 phagebacterial 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 $\approx $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 microcolonies or biofilms (HallStoodley 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 frequencydependent 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 $1H$ 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 ($1H$), lead to recurring epidemics. CRISPRbased 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 $\beta $ and latent period $\lambda $, and on the bacterial growth rate $\alpha $. 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 ${\delta}^{\star}$, 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 nonmicrobial systems. Our theoretical framework can be useful for identifying critical parameters, such as $H$ (and to some extent ${R}_{0}$), 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
Experimental methods
Engineering resistance
Request a detailed protocolOligonucleotides AAACTTCGGGAAGCACTTGTGGAAG and AAAACTTCCACAAGTGCTTCCCGAA were ordered from SigmaAldrich, 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 rfb50 rph1). The T7 wildtype phage was used in all experiments.
Efficiency of the CRISPRCas system
Request a detailed protocolEfficiency of the engineered CRISPRCas system was tested using the following protocol: Overnight cultures grown in LB containing 25 $\mu 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 $15\mathrm{min}$ in ${30}^{\circ}C$. Cells were spun down for 2 min in room temperature at $21130g$. Supernatant was discarded and the cell pellet was resuspended in 950 $\mu l$ of 1X TrisHCl buffer containing 0.4% ($\approx 227\mu M$) ascorbic acid prewarmed to ${43}^{\circ}C$ 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 $\mu g/ml$ chloramphenicol. These CRISPRCas system efficiencies at different MOIs were tested if they are statistically different from each other using twotailed unequal variances ttest at 0.05 confidence level using RStudio (R Core Team, 2013).
Determining the mean number of phages per cell
Request a detailed protocolThe cultures that were plated using standard plaque assays in the “Efficiency of the CRISPRCas system” experiment were also plated on LB agar plates containing 25 $\mu 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 protocolHerd immunity in a liquid culture was tested in 200 $\mu l$ of LB broth supplemented with 25 $\mu g/ml$ chloramphenicol in Nunclon flat bottom 96 well plate in a BioTek 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 $\approx {10}^{4}$ ($\approx 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.
Timelapse imaging of plaque growth
Request a detailed protocolSoft LB agar (0.7%) containing 25 $\mu g/ml$ chloramphenicol was melted and 3 $ml$ were poured into glass test tubes heated to ${43}^{\circ}C$ 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 $\mu l$ of bacteriophage stock, diluted such that there would be $\approx $10 plaques per plate, was added to the solution. Tubes were vortexed thoroughly and poured as an overlay on LB agar plates containing 25 $\mu 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 ${30}^{\circ}C$. 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. Timelapse 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 protocolGrowth 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 Accujet pro pipettor into 1 $ml$ of M9 buffer prewarmed to ${43}^{\circ}C$, vortexed for 15 seconds and incubated for 10 minutes in ${43}^{\circ}C$ with two more vortexing steps after 5 and 10 minutes of incubation. Samples were serially diluted and plated on LB agar plates containing 25 $\mu 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 protocolNutrientdependent growth rate of susceptible bacteria was measured in Nunclon flat bottom 96 well plate in BioTek Synergy H1 Plate reader in ${30}^{\circ}C$. 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 $\mu 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,
Results for the two fitting parameters, ${\alpha}_{\mathrm{max}}$ and ${K}_{c}$, 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 nutrientdependent growth rate measurements. Twosample ttest 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 $\mu g/ml$ chloramphenicol.
Phage burst sizes
Request a detailed protocolPhage burst sizes in bacteria growing at different growth rates were measured by onestep 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 ttest, $p\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0.05$). All experiments were performed in triplicates.
Phage latent periods
Request a detailed protocolPhage 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 protocolThe 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 protocolSoft M9 salts soft agar (0.5%) was supplemented with SYBR safe staining (final conc. 1%) and poured into glass cuvettes (VWR type 6040OG) to fill $\sim 2cm$ of the cuvette height. After soft agar solidification, the same stained soft agar was supplemented with T7 phage particles to a final concentration ${10}^{11}\mathrm{pfu}/\mathrm{ml}$ and poured on top of the agar without phages. The cuvettes were monitored in ${30}^{\circ}C$ 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 ${L}_{i}$ of fluorescence (a grayscale 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 ${L}_{i}$ 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 $\sim 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$,
Such a convolution with the heat kernel ${K}_{ij}(D)={(4\pi D)}^{1/2}\mathrm{exp}\left({(ij)}^{2}/4D\right)$ 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 ${\mathrm{pixel}}^{2}/\mathrm{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
which is in agreement with previous measures of phage diffusion (Stent and Wollman, 1952; Bayer and DeBlois, 1974; Briandet et al., 2008).
Modelling
Phage growth
Request a detailed protocolIn 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,
Both bacterial populations ${B}_{i}$, $i\in \{\mathrm{s},\mathrm{r}\}$, grow with rate $\alpha $ and decay via adsorption of phages $A[{B}_{i},P{B}_{\mathrm{s}},{B}_{\mathrm{r}}]$, an expression that is specified below. Infected populations ${I}_{i}$ gain numbers by adsorption and decrease via bursting. Resistant bacteria also can recover from their infected state with a recovery rate $\rho $. 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[{B}_{i},P{B}_{\mathrm{s}},{B}_{\mathrm{r}}]$, 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),
with an adsorption constant $\delta $. 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 LotkaVolterra 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 ${I}_{\mathrm{r}}$. 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 $\rho \gg 1/\lambda $. If furthermore $\rho \gg \delta 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, ${I}_{\mathrm{r}}\to 0$, 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 $\rho $.
Exponentially growing bacteria lead to double exponential phage growth
Request a detailed protocolFor convenience, we transform the populations to the total bacterial density $B={B}_{\mathrm{s}}+{B}_{\mathrm{r}}$ and introduce the fraction of susceptible cells $S={B}_{\mathrm{s}}/B$. The crucial assumption for the remainder of this section is that phages burst immediately after infection, $\lambda =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
If we assume that in initial stages of phage growth the number of phages is small, ie. $\delta P\ll \alpha \sim \mathcal{O}(1{h}^{1})$, the dynamics of bacteria and the fraction of susceptibles simplify to ${\partial}_{t}B=\alpha B$ and ${\partial}_{t}S=0$. Note that this term $\delta P$ also occurs in the linear phage dynamics, where it cannot be neglected. In this instance, we need to view $\delta B$ as a coefficient, which is likely much larger initially. This set of simplified equations can be solved in closed form,
The structure of phage dynamics is particularly important here – it exhibits a doubleexponential dependence on time $t$, which is a very fast, almost explosive, growth. Such doubleexponential 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
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 timedependence 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 doubleexponential growth of phages occurs. Hence, we compute the time ${T}_{\delta}$ defined as when phages reach a population of $P({T}_{\delta})=1/\delta $ assuming phages grow as Equation 15 until then. After ${T}_{\delta}$ the assumptions that allowed to obtain Equation 15 are not valid anymore. Inverting Equation 15 for time leads to
Subsequently, we can compare this estimate ${T}_{\delta}$ to the depletion time ${T}_{\mathrm{depl}}=(1/\alpha )\mathrm{log}({B}_{\mathrm{\infty}}/{B}_{0})$. When rearranging the inequality $T}_{\mathrm{d}\mathrm{e}\mathrm{p}\mathrm{l}}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{T}_{\delta$ in terms of the (initial) fraction of susceptibles ${S}_{0}$, we obtain
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}_{\mathrm{\infty}}$ usually fulfills $\delta {B}_{\mathrm{\infty}}/\alpha \gg 1$, such that the correction given by the second term of Equation 18 can be considered small. Thus, if phages grow ($\beta \phantom{\rule{thinmathspace}{0ex}}{S}_{0}\phantom{\rule{thinmathspace}{0ex}}>1$), they also grow very fast with a doubleexponential timedependence 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 protocolThe analysis above only treated the case $\lambda \to 0$. However, we reported that the latency time $\lambda $ increases significantly when bacterial growth rate $\alpha $ declines, see Table 2. Considering finite latency times entails dealing with an infected bacterial population $I$. (However, we identify $I\equiv {I}_{\mathrm{s}}$ and set ${I}_{\mathrm{r}}=0$.)
To this end, note that we can rearrange Equation 12 to $\left(1+\lambda {\partial}_{t}\right)I=\lambda \delta SBP$ using the adsorption model in Equation 13. Hence, we can use the differential operator $(1+\lambda {\partial}_{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
where we also inserted ${\partial}_{t}B\approx \alpha B$ in the last term, as we aim again for a solution at initial times where $\delta P\ll \alpha $. The effects of the limit $\lambda \to 0$ are directly observable – no terms are undefined in this limit. In particular, we find that equation Equation 19 and $\lambda =0$ lead directly to the dynamics of phages we just analyzed above, obtaining solution Equation 15.
In principle, Equation 19 is a hyperbolic reactiondiffusionequation, which is known to occur upon transformation (or approximation) of timedelayed differential equations (Fort and Méndez, 2002b). For initial times we can use the solutions $B(t)={B}_{0}\mathrm{exp}(\alpha t)$ and $S(t)={S}_{0}$. To proceed, we introduce the auxiliary variable
and assume $P(z)$ as a function of this new variable $z$. We need to transform the differential operators of time derivatives, and obtain ${\partial}_{t}=\frac{\partial z(t)}{\partial t}{\partial}_{z}=\alpha z{\partial}_{z}$ and ${\partial}_{t}^{2}=(\alpha z{\partial}_{z})(\alpha z{\partial}_{z})={\alpha}^{2}(z{\partial}_{z}+{z}^{2}{\partial}_{z}^{2})$. Inserting these expressions in Equation 19 and multiplying the whole equation with ${({\alpha}^{2}\lambda z)}^{1}$ yields the dynamics for phages,
where the two extant constants are $a=1(\beta {S}_{0}1)/(\lambda \alpha )$ and $b=1+1/(\lambda \alpha )$. Equation 21 is called 'Kummer’s equation' with confluent hypergeometric functions ${}_{1}F_{1}$ as solutions (Abramowitz and Stegun, 1964, pg. 504),
The two integration constants $A$ and $B$ can be determined via the initial conditions $P(t=0)={P}_{0}$ and $({\partial}_{t}P)(t=0)=\delta {B}_{0}{P}_{0}$. Using these conditions, the shape of the solution is again similar to before with $\lambda =0$ (double exponential timedependence), although $\lambda \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}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,
which is a nontrivial extension including finite latency times $\lambda $.
Note, however, that this relation Equation 23 does not indicate a correction to the general epidemiological parameter ${R}_{0}$, which can be identified with $\beta $ 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 $\lambda \alpha $ 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 $\alpha \to 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 $\alpha \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$. The most important of these is the transformation to $z(t)=\delta B(t)/\alpha $, which actually exhibits two problems: dividing by $\alpha $ 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 ${\partial}_{t}B=(\alpha \delta SP)B$ throughout our calculation. For $\alpha =0$ this second term is dominant in bacterial dynamics and would generate nonlinear phage dynamics if inserted for ${\partial}_{t}B$ 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 $\beta {S}_{0}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1$ could indicate phage expansion and bacterial decay.
Growth of phages on plated bacterial lawn
Request a detailed protocolSpatial 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; OrtegaCejas 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, ${B}_{i}={B}_{i}(\overrightarrow{x},t)$, $i\in \{\mathrm{s},\mathrm{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 $\widehat{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, $\widehat{B}(B)\approx const$. Consequently, the adsorption term can be written in the following way,
which only depends on the relative frequencies of bacterial strains. The adsorption constant ${\delta}^{\star}$ is both the rate of adsorption and interhost transit time as determined by the diffusion constant $D$. Thus, one can expect the formal dependence ${\delta}^{\star}={\delta}^{\star}(D,\widehat{B}(B))$. For our particular experimental setup, however, ${\delta}^{\star}$ will be treated as a constant. This adsorption term Equation 24 leads to the dynamics of phages
where we collected all contributions to phage growth in $G[P,S]$ and added the spatial diffusion term $D{\nabla}^{2}P$. For simplicity, we consider only expansion in a single dimension (${\nabla}^{2}\equiv {\partial}_{x}^{2}$), 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,
where we also consider the correction $\lambda \alpha $ 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. Reactiondiffusion 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,
Only the linearized growth rate of phages at very low densities is relevant for the expansion speed, ${\partial}_{P}G[P=0,S]$. Thus, the fraction of susceptible individuals $S$ should be unchanged from its initial value ${S}_{0}$. It should be noted, that only for ${S}_{0}\beta \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1+\lambda \alpha$ 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 ${\beta}_{\mathrm{depl}}$ and ${\lambda}_{\mathrm{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 CRISPRCas 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 unstructured 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 ${T}_{\mathrm{depl}}$ 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)={\int}_{0}^{t}d{t}^{\prime}v({t}^{\prime})$. Employing the simplification that only two values of phage growth are necessary to describe the dynamics – before ${T}_{\mathrm{depl}}$ phages grow normally with $\beta $ and $\lambda $, after ${T}_{\mathrm{depl}}$ phage growth changes to ${\beta}_{\mathrm{depl}}$ and ${\lambda}_{\mathrm{depl}}$ – we can evaluate the integral for the radius directly, arriving at,
Using this expression we estimated the adsorption constant ${\delta}^{\star}$ from the growth experiments as it difficult to measure in practice. This estimate is done for radii exactly at the time of nutrient depletion ${T}_{\mathrm{depl}}$, 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 ${\beta}_{\mathrm{depl}}=3.0$, which results in ${H}_{\mathrm{depl}}=\left({\beta}_{\mathrm{depl}}1\right)/{\beta}_{\mathrm{depl}}\approx 0.67$. Thus, all experiments with $S\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}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 ($S\le 0.9$), which would correspond to ${\beta}_{\mathrm{d}\mathrm{e}\mathrm{p}\mathrm{l}}\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}1.1$. This value is, however, still within experimental accuracy of our estimates of ${\beta}_{\mathrm{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 $\rho $ (instead of the $\rho \to \mathrm{\infty}$ 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 $\delta $, the recovery rate $\rho $ and the yield coefficient $Y$. Undiluted LB medium is known to support a population of $5\cdot {10}^{9}\mathrm{cells}/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 $\delta $ 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 $\delta $ is changed by orders of magnitude, $\delta \sim \mathcal{O}\left({10}^{6}\mathrm{\dots}{10}^{8}\right)$. 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 $\delta ={10}^{7}{h}^{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 $\rho $. A first indication of the value of $\rho $ 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/\lambda \sim \mathcal{O}(1)$ to estimate $\rho \sim \mathcal{O}\left({10}^{3}\right)$. However, our results also indicate that recovery via the CRISPR/Cas system heavily depends on MOI, implying that $\rho $ 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 $\rho $ 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 $\rho \sim \mathcal{O}\left(1\right)$. Lower values of $\rho $ do not allow the resistant population to recover from phage infection, while for larger values of $\rho $, phages are drained from the culture very fast. Such a small value of $\rho $ 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 $\rho =1.5{h}^{1}$ is the recovery rate when the CRISPR/Cas system is heavily stressed, which is comparable to the actual burst rate $1/\lambda $ 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, nonmonotonic 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.
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 $P\sim 1/\delta $, 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 $P\sim 1/\delta $ is exceeded.
In order to proceed, we investigate the system at time ${T}_{\delta}$ 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\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1$ infections before they burst and produce only $\kappa \beta $ phages, which reduces the burst size by a (yet unspecified) factor $0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}\kappa <1$. $\kappa =1$ implies that resistant cells produce the same number of phages as susceptible cells, while $\kappa =0$ indicates only cell death. Combining these considerations yields
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}_{\delta})$ can be estimated by inserting the time ${T}_{\delta}$ from Equation 18 into the exponential growth Equation 15b. Subsequently, we can rearrange Equation 29, obtaining a bound on $M$:
The first term $1/\delta B({T}_{\delta})$ indicates the ratio of phages to bacteria at time ${T}_{\delta}$, and can be considered small for nonextremal parameters compared to the other terms. This fact justifies our assumption that the actual value of $\delta $ 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 $M\lesssim 3+86\kappa$. Thus, each resistant bacterial cell could degrade up to $\mathcal{O}\left({10}^{1}\mathrm{\dots}{10}^{2}\right)$ 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—source data 1
 https://doi.org/10.7554/eLife.32035.027
Appendix 3
Additional information on experimental setup
Our experimental setup for the scanner system is shown in Appendix 3—figure 1.
Appendix 4
Source data and code
Timelapse 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/elifesciencespublications/phagegrowth.
Data availability

Data from: CRISPRbased herd immunity limits phage epidemics in bacterial populationsAvailable at Dryad Digital Repository under a CC0 Public Domain Dedication.
References

Optimizing bacteriophage plaque fecundityJournal of Theoretical Biology 249:582–592.https://doi.org/10.1016/j.jtbi.2007.08.006

BookHandbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, 55Courier Corporation.

Spatial, temporal, and genetic heterogeneity in host populations and the design of immunization programmesMathematical Medicine and Biology 1:233–266.https://doi.org/10.1093/imammb/1.3.233

Diffusion constant and dimension of bacteriophage phi X174 as determined by selfbeat laser light spectroscopy and electron microscopyJournal of Virology 14:975–980.

Modern Bicrobial Genetics177–225, Restrictionmodification systems, Modern Bicrobial Genetics, 10.1002/047122197X.ch7.

Fluorescence correlation spectroscopy to study diffusion and reaction of bacteriophages inside biofilmsApplied and Environmental Microbiology 74:2135–2143.https://doi.org/10.1128/AEM.0230407

Immigration of susceptible hosts triggers the evolution of alternative parasite defence strategiesProceedings of the Royal Society B: Biological Sciences 283:20160721.https://doi.org/10.1098/rspb.2016.0721

Migration of chemotactic bacteria in soft agar: role of gel concentrationBiophysical Journal 101:525–534.https://doi.org/10.1016/j.bpj.2011.06.023

Smallpox: emergence, global spread, and eradicationHistory and Philosophy of the Life Sciences 15:397–420.

Network frailty and the geometry of herd immunityProceedings of the Royal Society B: Biological Sciences 273:2743–2748.https://doi.org/10.1098/rspb.2006.3636

"Herd immunity": a rough guideClinical Infectious Diseases 52:911–916.https://doi.org/10.1093/cid/cir007

Herd immunity: history, theory, practiceEpidemiologic Reviews 15:265–302.https://doi.org/10.1093/oxfordjournals.epirev.a036121

The wave of advance of advantageous genesAnnals of Eugenics 7:355–369.https://doi.org/10.1111/j.14691809.1937.tb02153.x

Timedelayed spread of viruses in growing plaquesPhysical Review Letters 89:178101.https://doi.org/10.1103/PhysRevLett.89.178101

Wavefronts in timedelayed reactiondiffusion systems. Theory and comparison to experimentReports on Progress in Physics 65:895–954.https://doi.org/10.1088/00344885/65/6/201

Herd immunity: basic concept and relevance to public health immunization practicesAmerican Journal of Epidemiology 94:179–189.https://doi.org/10.1093/oxfordjournals.aje.a121310

Rate of translocation of bacteriophage T7 DNA across the membranes of Escherichia coliJournal of Bacteriology 177:4066–4076.https://doi.org/10.1128/jb.177.14.40664076.1995

Demography and diffusion in epidemics: malaria and black death spreadActa Biotheoretica 58:277–305.https://doi.org/10.1007/s104410109103z

Mathematical models of infectious disease transmissionNature Reviews Microbiology 6:477–487.https://doi.org/10.1038/nrmicro1845

Bacterial biofilms: from the natural environment to infectious diseasesNature Reviews Microbiology 2:95–108.https://doi.org/10.1038/nrmicro821

BookEpidemic Disease in England: The Evidence of Variability and of Persistency of TypeBedford Press.

A brief history of R0 and a recipe for its calculationActa Biotheoretica 50:189–204.https://doi.org/10.1023/A:1016599411804

BookEvolutionary Games and Population DynamicsCambridge, UK: Cambridge University Press.https://doi.org/10.1017/CBO9781139173179

Pattern formation triggered by rare events: lessons from the spread of rabiesProceedings of the Royal Society B: Biological Sciences 264:495–503.https://doi.org/10.1098/rspb.1997.0071

RNAguided editing of bacterial genomes using CRISPRCas systemsNature Biotechnology 31:233–239.https://doi.org/10.1038/nbt.2508

BookA Study of the Diffusion Equation with Growth of the Quantity of Matter and Its Application to a Biology ProblemIn: Tikhomirov VM, editors. Selected Works of A. N. Kolmogorov, 1. Dordrecht, The Netherlands: Kluwer Academic Publishers. pp. 242–270.

Spatial heterogeneity in epidemic modelsJournal of Theoretical Biology 179:1–11.https://doi.org/10.1006/jtbi.1996.0042

Population dynamics of phage and bacteria in spatially structured habitats using phage λ and Escherichia coliJournal of Bacteriology 198:1783–1793.https://doi.org/10.1128/JB.0096515

Mechanism of Inactivation of Bacteriophage J1 by Ascorbic AcidAgricultural and Biological Chemistry 37:1145–1151.https://doi.org/10.1080/00021369.1973.10860809

On the spatial spread of rabies among foxesProceedings of the Royal Society B: Biological Sciences 229:111–150.https://doi.org/10.1098/rspb.1986.0078

BookMathematical Biology I: An Introduction, Vol. 17 of Interdisciplinary Applied MathematicsNew York: Springer.

Effect of protein A on adsorption of bacteriophages to Staphylococcus aureusJournal of Virology 14:198–202.

Approximate solution to the speed of spreading virusesPhysical Review E 69:031909.https://doi.org/10.1103/PhysRevE.69.031909

phagegrowthphagegrowth, 42ec97e, https://github.com/lukasgeyrhofer/phagegrowth.

A model of the velocity of advance of foot and mouth disease in feral pigsThe Journal of Applied Ecology 27:635–650.https://doi.org/10.2307/2404308

SoftwareR: a language and environment for statistical computingR Foundation for Statistical Computing, Vienna, Austria.

Spatial dynamics and genetics of infectious diseases on heterogeneous landscapesJournal of The Royal Society Interface 4:935–948.https://doi.org/10.1098/rsif.2007.1041

An agestructured model of pre and postvaccination measles transmissionMathematical Medicine and Biology 1:169–191.https://doi.org/10.1093/imammb/1.2.169

On the velocity of epidemic waves in model plant disease epidemicsEcological Modelling 87:217–222.https://doi.org/10.1016/03043800(95)000305

Fiji: an opensource platform for biologicalimage analysisNature Methods 9:676–682.https://doi.org/10.1038/nmeth.2019

CRISPRmediated adaptive immune systems in bacteria and archaeaAnnual Review of Biochemistry 82:237–266.https://doi.org/10.1146/annurevbiochem072911172315

On the two step nature of bacteriophage absorptionBiochimica et Biophysica Acta 8:260–269.https://doi.org/10.1016/00063002(52)900413

The interaction of phage and biofilmsFEMS Microbiology Letters 232:1–6.https://doi.org/10.1016/S03781097(04)000412

Rapidly evolving CRISPRs implicated in acquired resistance of microorganisms to virusesEnvironmental Microbiology 10:200–207.https://doi.org/10.1111/j.14622920.2007.01444.x

Front propagation into unstable statesPhysics Reports 386:29–222.https://doi.org/10.1016/j.physrep.2003.08.001

BookQuantitative Viral Ecology: Dynamics of Viruses and Their Microbial HostsPrinceton, United States: Princeton University Press.

A genetic herdimmunity model for the maintenance of MHC polymorphismImmunological Reviews 143:263–292.https://doi.org/10.1111/j.1600065X.1995.tb00679.x

Replication of viruses in a growing plaque: a reactiondiffusion modelBiophysical Journal 61:1540–1549.https://doi.org/10.1016/S00063495(92)819586

Seasonality and the requirements for perpetuation and eradication of viruses in populationsAmerican Journal of Epidemiology 109:103–123.https://doi.org/10.1093/oxfordjournals.aje.a112666

Amplification and spread of viruses in a growing plaqueJournal of Theoretical Biology 200:365–373.https://doi.org/10.1006/jtbi.1999.1001
Decision letter

Richard A NeherReviewing 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 "CRISPRbased 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.
Summary:
This manuscript presents a mathematical model and a set of experiments examining how CRISPRCas mediated immunity to phage infection can result in herd immunity preventing the spread of the pathogen. The experiments are performed both in wellmixed 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 longrange jumps ultimately to the wellmixed 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 crossover 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.
https://doi.org/10.7554/eLife.32035.034Author 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 prerequisite 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 longrange jumps ultimately to the wellmixed 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 longrange 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 longrange 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.
https://doi.org/10.7554/eLife.32035.035The 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 crossover 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
Funding
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.
Acknowledgements
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
 Richard A Neher, University of Basel, Switzerland
Publication history
 Received: October 7, 2017
 Accepted: March 8, 2018
 Accepted Manuscript published: March 9, 2018 (version 1)
 Version of Record published: April 27, 2018 (version 2)
 Version of Record updated: June 11, 2018 (version 3)
Copyright
© 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.
Metrics

 2,827
 Page views

 424
 Downloads

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