CRISPR-based herd immunity limits phage epidemics in bacterial populations

Pathogens are a driving force in evolution, giving rise to a diversity of host immune defenses. In order for a pathogen to spread in a population a sufficient number of its members must be susceptible to infection, as resistant individuals can prevent the spread of a pathogen among susceptible hosts in a process known as herd immunity. While herd immunity has been extensively studied in vertebrate populations, little is known about its role, if any, in the dynamics between bacteria and their phage pathogens. Here we explore the dynamics of T7 phage epidemics in structured and unstructured Escherichia coli populations consisting of differing mixtures of susceptible and resistant individuals harboring CRISPR immunity to the phage. Using both experiments and mathematical modelling we describe the conditions under which herd immunity arises in bacterial populations. Notably, the effects of herd immunity depend strongly on the presence of spatial structure in the population, the bacterial growth rate, and phage replication rate. The results of our model can apply to other host-pathogen systems to determine the herd immunity threshold from the relative speed of an epidemic wave in partially resistant populations. In addition, our findings suggest that herd immunity plays an important role in bacterial communities, as seen in other host-pathogen systems, allowing for stable coexistence of bacteria and their phages and the maintenance of polymorphism in bacterial immunity.

Cas9 Pre-crRNA transcript  Figure 1. Mechanism of CRISPR/Cas type II immunity. The CRISPR/Cas system provides immunity to phages and its main features can be described by three distinct stages. (A) Acquisition. When a cell gets infected by a phage, a protospacer on the invading phage DNA (indicated as a red bar) is recognized by Cas1 and Cas2.
The protospacer is cleaved out and ligated to the leader end (proximal to the Cas genes) of the CRISPR array as a newly acquired spacer (red diamond). (B) Processing. The CRISPR array is transcribed as a Pre-crRNA and processed by Cas9 (assisted by RNaseIII and trans-activating RNA, not shown) into mature crRNAs. (C) Interference. Mature crRNAs associate with Cas9 proteins to form interference complexes which are guided by sequence complementarity between the crRNAs and protospacers to cleave invading DNA of phages whose protospacers have been previously incorporated into the CRISPR array. (D) A truncated version of the CRISPR system on a low copy plasmid, which was used in this study lacks cas1 and cas2 genes and was engineered to target a protospacer on the T7 phage chromosome to provide Escherichia coli cells with immunity to the phage. The susceptible strain contains the same plasmid except the spacer does not target the T7 phage chromosome.
is defined as, which determines the required minimum fraction of resistant individuals needed to halt the spread 42 of an epidemic and is effected by the specific details of transmission and the contact rate among 43 individuals Grassly and Fraser (2008). Many theoretical studies have addressed the influence of 44 some of these details, in particular maternal immunity Anderson and May (1992), age at vaccination 45 Anderson and May (1982); Nokes and Anderson (1988), age related or seasonal differences in 46 contact rates Schenzle (1984); Anderson and May (1985); Yorke et al. (1979), social structure Fox 47 are required. 79 To investigate under which conditions herd immunity may arise in bacterial populations, we 80 constructed an experimental system consisting of T7 phage and bacterial strains susceptible and 81 resistant to it. Our experimental system can be characterized by the following features. First, 82 we used two strains of Escherichia coli, one with an engineered CRISPR-based immunity to the T7 83 phage, and the other lacking it ( Fig 1D). Second, we examined the dynamics of the phage spread in 84 different environments -spatially structured and without structure. Furthermore, we developed 85 and analyzed a spatially explicit model of our experimental system to determine the biologically 86 relevant parameters necessary for bacterial populations to exhibit herd immunity. 87

88
Properties of resistant individuals 89 We engineered a resistant E. coli strain by introducing the CRISPR-Cas Type II system from Streptococ-90 cus pyogenes with a spacer targeting the T7 phage genome (see Material and Methods). We further 91 characterized the ability of the system to confer resistance to the phage. We find a significant level 92 of resistance as measured by the probability of cell burst when exposed to T7 (Fig 2A). However, 93 resistance is not fully penetrant as approximately 1 in 1000 resistant cells succumb to infection. In 94 addition, we observe that as phage load increases (multiplicity of infection, MOI) the probability 95 that a cell bursts increases (Fig 2A). In order to determine the herd immunity threshold in our 96 experimental system, we constructed the resistant strain such that upon infection the cell growth 97 is halted, yet the cell still adsorbs and degrades phages (Fig 2B,C). This feature is important as it 98 prevents the action of frequency dependent selection which in naturally growing populations will 99 favor the resistant strain until its frequency reaches the herd immunity threshold. Thus, in our 100 system if the frequency of the resistant strain is below the herd immunity threshold, the resistant 101 cells remain below the threshold and are unable to stop the epidemic and the whole population 102 collapses. In contrast, if the frequency of resistant individuals in the population is above the herd 103 immunity threshold, the resistant individuals provide complete herd immunity and the population 104 survives. These properties allow us to quantify the expanding epidemic in both liquid media and on 105 bacterial lawns (without and with spatial structure, respectively) using high throughput techniques.    Table 2 for burst size estimates). Each population phage challenge is replicated 16 times. The solid dark green line shows the model prediction, Eqn. (4), for the herd immunity threshold ( ), given latent period ( ), bacterial growth rate ( ), and phage burst size ( ). Shaded area indicates ±1 standard deviation.
Specifically, it allows us to control for the complex dynamics of the system arising from frequency 107 dependent selection and simultaneous changes in the physiological states of the cells (growth rates 108 depending on the nutrient concentrations) and phage (burst size, latent period depending on the  Table 2. (C) Population survival analysis upon phage challenge as a function of the fraction of resistant cells and the intrinsic growth rate (nutrient availability, ). The empty circles, circled dots, and filled circles represent the outcomes of 18 replicates done in 3 independent batches, as indicated in the legend. The green line shows the model predicted herd immunity threshold given by Eqn. (5), using a quadratic fit for ∕ and inverting the Monod kinetics of bacterial growth (see Fig 6B)  and collapse under the epidemic (Fig 3). 125 As mentioned in the introduction, phage and bacterial physiology may affect the herd immunity 126 threshold. To test this we altered bacterial growth by reducing the concentration of nutrients in the 127 medium (Fig 6) which concurrently alters the T7 phage's latent period and burst size (Fig 4A,   susceptible individuals can be easily quantified. In addition, these data allow for estimating the 148 speed of the epidemic wave front in these different regimes using real-time imaging (Fig 5A). 149 We observe a decline in the number of plaque forming units (data not shown) and a significant the plaques are visually detectable (Fig 7). It should be noted that all populations with such low 162 percentages of resistant individuals in liquid environment collapsed, indicating that spatial structure 163 significantly facilitates herd immunity. 164 The results presented in this and the previous section would allow us to use Eq. (1) to infer a 165 value for 0 from the observed threshold between surviving and collapsing bacterial populations, 166 Figs 3 and 4. We also observe that herd immunity is strongly influenced by spatial organization of 167 the population, Fig 5. How the exact value of (and subsequently the "classical" epidemiological 168 parameter 0 ) is affected by various factors such as bacterial growth rate, phage burst size and 169 latent period is, however, difficult to resolve experimentally. Similarly, quantification of the effect of 170 spatial structure is hardly achievable solely by experimental investigation. In order to disentangle 171 the roles of all the factors mentioned above, we proceed with development and analysis of a 172 mathematical model of the experimental system.

173
Modelling bacterial herd immunity 174 We developed a model of phage growth that takes several physiological processes into account: 175 bacterial growth during the experiment, bacterial mortality due to phage infection, and phage 176 mortality due to bacterial immunity. Furthermore, we use the previously reported observation that 177 phage burst size and latent period depend strongly on the bacterial growth rate (see Table 1). 178 The main processes in our model system can be defined by the following set of reactions, and burst size decreases, such that phage reproduction declines dramatically (see Table 2). We ( Here, 0 and ∞ are the initial and final bacterial densities, respectively. In the initial exponential 193 growth phase, our estimates from experimental data for growth parameters are = 0.63 ℎ −1 , = 194 85.6 phages∕cell and = 0.60 ℎ, for bacteria and phages, respectively (Table 1 and Table 2 Table 2). From 225 these data we estimated the dependence of the phage burst size on the bacterial growth rate, ( ), 226 using a numerical quadratic fit ( Fig 4A). Similarly, we estimated the dependence of the phage latent 227 period on the bacterial growth rate, ( ) ( Fig 4B). Using these estimates we calculated the expected 228 growth rate-dependent herd immunity threshold which gives a very good prediction of the shift in the herd immunity threshold to lower values 230 for slower growing populations (green line in Fig 4C). This verification of our model shows that it 231 correctly captures the dependence of the herd immunity threshold on bacterial and phage growth 232 parameters.

233
Modelling herd immunity in spatially structured populations 234 The dynamics of phage spread differ between growth in unstructured (e.g., liquid) and structured 235 (e.g., plates) populations. In order to quantify the effect of spatial structure in a population, we 236 extend our model to include a spatial dimension. In structured populations growth is a radial which is computed in more details in the Materials and Methods. This expression (6) indicates that 248 the population composition crucially influences the spreading speed at much lower fractions of 249 resistant bacteria than the herd immunity threshold (4), where phage expansion stops completely. 250 Consequently, the resulting plaque radius decays with increasing fractions of resistants and 251 reaches zero at . A prediction for can be obtained by integrating (6) This expression, (7), is devoid of any (explicit) environmental conditions, which are not already 260 contained in the herd immunity threshold itself. Thus, it could apply to any pathogen-host shown to facilitate herd immunity. These suggest that herd immunity may be fairly prevalent in low 286 nutrient communities such as soil and oligotrophic marine environments. 287 In an evolutionary context, herd immunity may also impact the efficacy of selection as the 288 selective advantage of a resistance allele will diminish as the frequency of the resistant allele in 289 a population approaches the herd immunity threshold, . This has two important implications. 290 First, while complete selective sweeps result in the reduction of genetic diversity at linked loci, herd 291 immunity may lead to only partial selective sweeps in which some diversity is maintained. Second, 292 herd immunity has a potential to generate and maintain polymorphism at immunity loci, as has been 293 shown for genes coding for the major histocompatibility complex (MHC) Wills and Green (1995).   302 We also developed a mathematical model and show how the herd immunity threshold 303 (Eqn. (4)) depends on the phage burst size and latent period , and on the bacterial growth rate . 304 This dependence arises as phages have to outgrow the growing bacterial population, as host and 305 pathogen have similar generation times in our microbial setting. In addition to these parameters, 306 we also describe how the speed (Eqn. (6)     LB 0 0.0 101.1 (± 10.9) 3.0 (± 1.9) 1.8 (± 1.1) LB 20 0.2 43.4 (± 3.9) 12.0 (± 4.2) 16.6 (± 6.0) LB 50 0.5 40.0 (± 3.0) 35.6 (± 16.4) 53.4 (± 24.9) LB 100 1.0 36.1 (± 6.1) 85.6 (± 47.3) 142.1 (± 82.1) If we assume that in initial stages of phage growth the number of phages is small, ie. ≪ ∼ 486 (1 ℎ −1 ), the dynamics of bacteria and the fraction of susceptibles simplify to = and = 0.
Note that this term also occurs in the linear phage dynamics, where it cannot be neglected. In 488 this instance, we need to view as a coefficient, which is likely much larger initially. This set of 489 simplified equations can be solved in closed form, The structure of phage dynamics is particularly important here -it exhibits a double-exponential This expression (17) is a condition for phages to reach "large" population sizes before nutrients 508 are depleted by bacteria. The final population density ∞ usually fulfills ∞ ≫ 1, such that the 509 correction given by the second term of (17) can be considered small. Thus, if phages grow ( 0 > 1), 510 they also grow very fast with a double-exponential time-dependence and reach a considerably large 511 population size before bacteria stop multiplying (for almost all parameter values).

512
Extending analysis to finite burst times 513 The analysis above only treated the case → 0. However, we reported that the latency time 514 increases significantly when bacterial growth rate declines, see Table 2. Considering finite latency 515 times entails dealing with an infected bacterial population . (However, we identify ≡ s and set 516 r = 0.)

517
To this end, note that we can rearrange (11a) to 1 + = using the adsorption model 518 in (12). Hence, we can use the differential operator (1+ ) and apply it directly to (11e) to reduce the 519 dependence on in this equation at the cost of introducing higher order derivatives. In particular, 520 we obtain where we also inserted ≈ in the last term, as we aim again for a solution at initial times 522 where ≪ . The effects of the limit → 0 are directly observable -no terms are undefined in this limit. In particular, we find that equation (18) and = 0 lead directly to the dynamics of phages 524 we just analyzed above, obtaining solution (14c). 525 In principle, (18) is a hyperbolic reaction-diffusion-equation, which is known to occur upon 526 transformation (or approximation) of time-delayed differential equations Fort and Méndez (2002b). 527 For initial times we can use the solutions ( ) = 0 exp( ) and ( ) = 0 . To proceed, we introduce 528 the auxiliary variable and assume ( ) as a function of this new variable . We need to transform the differential operators which is a non-trivial extension including finite latency times .
where we collected all contributions to phage growth in [ , ] and added the spatial diffusion 582 term ∇ 2 . For simplicity, we consider only expansion in a single dimension (∇ 2 ≡ 2 ), which has 583 been found to coincide well with the dynamics of plaque growth Yin and McCaskill (1992). The 584 growth term for phages is then defined as, where we also consider the correction obtained from the analysis in liquid culture. Reaction- Only when nutrients are depleted, can pure diffusion processes explain the slow decrease in speed 607 observed in experiments (see Fig 7A). Our model assumes a sharp drop to = 0 at depl for small . 608 In order to derive an expression for the plaque radius , we integrate the expansion speed (6) 609 over time, ( ) = ∫ 0 ′ ( ′ ). Employing the simplification that only two values of phage growth are 610 necessary to describe the dynamics -before depl phages grow normally with and , after depl 611 phage growth changes to depl and depl -we can evaluate the integral for the radius directly, arriving 612 at, assume that is constant. 838 We ran parameter sweeps in simulations and compared the outcome -collapsed or surviving 839 bacterial populations -to the observed experimental results (see Fig 3). The best agreement  Table 3  to produce enough phages and ultimately the whole bacterial population goes extinct. 855 The purpose of the extended model in this section was to justify the fact that phages can wipe 856 out the whole bacterial population, which was not possible in the simplified model used in the main 857 text. There, the resistant bacterial population was basically unaffected by phages and just acted as Infection load and efficiency of the CRISPR/Cas system 861 In the section Modelling we showed that positive phage growth leads eventually to a very fast 862 increase in the phage population, that occurs before nutrients are depleted (for almost all realistic 863 parameters). This behavior of the dynamics was also observed in the extended simulation model 864 presented in the last section. Moreover, as a condition we used that the phage population reaches 865 a size ∼ 1∕ , which is after all arbitrary -it only determines if we can employ useful simplifications 866 and approximations to model equations. However, simulation results presented in the last section 867 indicate that the bacterial population starts to decay soon after such a threshold ∼ 1∕ is 868 exceeded. 869 In order to proceed, we investigate the system at time further. We assume that the phage 870 population is large enough that it will not be degraded by the CRISPR/Cas immune system. The 871 threat to immediate phage extinction is low at this point. The actual equations are hard to solve 872 directly, hence we revert to simple balance equations, ignoring the dynamical component. Specif- 873 ically, we compare the number of (present and eventually produced) phages to the number of 874 infections needed to wipe out the whole population. To incorporate the effects of the bacterial 875 immune system in resistant bacteria, we assume that they need > 1 infections before they burst 876 and produce only phages, which reduces the burst size by a (yet unspecified) factor 0 < < 1. Appendix: Additional information on experimental setup 892 Our experimental setup for the scanner system is shown in Fig. 10. 893