Ecological memory preserves phage resistance mechanisms in bacteria

Bacterial defenses against phage, which include CRISPR-mediated immunity and other mechanisms, can carry substantial growth rate costs and can be rapidly lost when pathogens are eliminated. How bacteria preserve their molecular defenses despite their costs, in the face of variable pathogen levels and inter-strain competition, remains a major unsolved problem in evolutionary biology. Here, we present a multilevel model that incorporates biophysics of molecular binding, host-pathogen population dynamics, and ecological dynamics across a large number of independent territories. Using techniques of game theory and non-linear dynamical systems, we show that by maintaining a non-zero failure rate of defenses, hosts sustain sufficient levels of pathogen within an ecology to select against loss of the defense. This resistance switching strategy is evolutionarily stable, and provides a powerful evolutionary mechanism that maintains host-pathogen interactions, selects against cheater strains that avoid the costs of immunity, and enables co-evolutionary dynamics in a wide range of systems.


Supplementary Figures
Supplementary Figure 1. A legend showing how to read flow diagrams as triangle plots. Corners correspond to pure populations of sensitive bacteria (S), resistant bacteria (R), or phage (P ). Frequencies of resistant, sensitive, and infected host and phage increase in the direction of the arrows along the edges of the simplex. Interior points correspond to mixtures of bacteria and phage, whose frequencies can be read by following dotted lines as illustrated in this example.
Supplementary Figure 2. Plot showing the growth rate of host cells and the phenotype and phage frequencies in the model of preventative defenses as α is varied across γ and γ at burst size β = 10. The system goes through a series of two transcritical bifurcations. S, SP and RSP denote stable phases of phage-host dynamics, corresponding to the fixation of sensitive bacteria (S), coexistence of sensitive bacteria and phage (SP), or coexistence of resistant and sensitive host and phage (RSP). Solid curves -stable fixed points, dashed curves -unstable fixed points. Remaining parameters are nr = 1, d = 1, b = 0.9, kL = 1, s = 0.01, Km = 0. From Eq. 12 we observe that when the RSP fixed point is stable, the growth rate of host cells is λ = b − s. For stable fixed points with fR = 0, we useḟS = 0 to obtain λ = d in the S phase and λ = (d−αnr)(βkL −d−kL)/(βkL −(1+nr)(d+kL)) in the SP phase. For periodic dynamics in the RSP phase, at steady state the average growth rate can be obtained by integratinġ fR/fR over the period T : 1 Since fR(t + T ) = fR(t), we obtain λ T = b − s, the growth rate of the resistant phenotype.
Supplementary Figure 3. Phase diagrams for preventative defenses for larger switching rates (panels a and b) and larger costs of defense mechanism, d − b (panel b). Insets show the separation of γ and γ for α > d. The RSP&S bistable region located between γ and γ curves becomes larger with the increase in switching rate and cost d − b. Remaining parameters are nr = 1, d = 1, kL = 1, Km = 0.
Supplementary Figure 4. Phase diagrams without resistance, as a function of the number of receptors nr, phage burst rate βkL and phage absorption rate α. Panels (a-d) correspond to increasing numbers of receptors on the cells, from nr = 1 (panel a), nr = 100 (panel b), nr = 1, 000 (panel c) to the limiting value of nr = ∞ (panel d). The distinct stable outcomes correspond to the S, SP and E phases, as defined in the text. Vertical dotted line indicates α = d/nr. Host extinction is stable for α > d/nr. Parameters used are d = 1, kL = 1 and Km = 0.
Supplementary Figure 5. Phase diagrams for two phage types, P1 and P2, shown for P2 parameters in the model of preventative defenses. (a-c) P1 parameters (α1, β1) are located in the RSP phase. Dashed curves separate regions of single-and multi-phage phases. The pink region corresponds to the case where P2 drives P1 and the resistant phenotype extinct. The resistant phenotype switches at a rate s = 10 −3 . Size of the coexistence regions weakly depends on s. Remaining parameters are d = 1, b = 0.9, kL = 1, Km = 0. We consider the phase structure of Eq. (12) for a range of nr including nr 1. In the limit nr → ∞ the dependence on nr drops out, as limn r→∞ kI = αP/A. As nr increases the SP phase moves to lower values of α < d/nr while the RSP phase for immune defenses develops a larger region where periodic dynamics are possible (shown in white). In the limit nr → ∞ the SP phase disappears and the phase diagram contains only the S and RSP phases separated by γ = (d + α)(d + kL)/α. Parameters used to generate the figure are d = 1, b = 0.9, s = 10 −3 , Km = 0. Here we consider phage that does not fall off the receptor upon infection and therefore blocks the receptor for subsequent absorption of phage in the infected phenotype. This effect would be small for large nr and most extreme in the model of minimal sensitivity (nr = 1), which would prevent the infected phenotype to act as a phage sink. We solved the nr = 1 model to find this change shifts the location of γ and γ bifurcations as well as impacts the dynamics at the intersection of S, SP and RSP phases. Parameters used are nr = 1, d = 1, b = 0.9, kL = 1, s = 10 −3 , κ = 0.

A. Turbidostat control of host biomass
We consider the control of host biomass by a feedback dilution mechanism that maintains a fixed host concentration within the system volume. The dynamics are given by transforming the simplex equations (12) to a set of coordinates given by phenotype frequencies in the host population x = R/(R+S+I), y = S/(R+S+I), z = I/(R+S+I) = 1−x−y and phage pressure p = P/(R + S + I). In this basis, the infection rate is k I = αn r p/(K m + n r a + p), with K m ≡ K m /(R + S + I) and where a = A/(R + S + I) is total frequency of phage-adsorbing hosts, which equals 1 − x for preventative and 1 for immune defenses. The turbidostat equations are: where λ(t) = dy + bx − k L z is the instantaneous host growth rate. These coordinates can be mapped to the simplex via: except at the host extinction point E, which is a singular outcome in this basis and needs to be studied on the simplex. The stability of fixed points in the host biomass basis matches that on the simplex and generates the same phases as the ones presented for the simplex equations.

B. Chemostat control of host biomass
We analyze the dynamics in a chemostat model given by the equationṡ Here ρ is the resource, η is the resource to growth conversion factor (which we set to 1), ρ in is the concentration of the supplied resource while D is the chemostat dilution rate. R, S, I and P denote concentrations of host and phage, and the infection rate is k I = αn r P/(K m + n r A + P ). Initially, we consider the high binding affinity limit and set K m = 0. Cellular division rates depend on resource concentration and are generally chosen to be Monod functions of the form d(ρ) = d ρ/(ρ c + ρ), where ρ c is nutrient concentration at half-maximal growth rate. We consider the regime of nutrient limitation ρ ρ c , for which division rates linearly depend on ρ: d(ρ) = d · ρ/ρ c , and similarly b(ρ) = b · ρ/ρ c ; we absorb ρ c into the prefactor by setting ρ c = 1. Extinction of both phage and host is globally stable if D > d · ρ in , which means the population is washed out faster than it can grow. Dilution rates below d · ρ in but above b · r in will select against the resistant phenotype. We will consider dilution rates D < b · ρ in − s so that the chemostat supports growth of all phenotypes.
The analysis of chemostat fixed points gives the same phases as in the turbidostat formulation with: (i) the E phase, which is unstable for our choice of dilution rate; (ii) the S phase, corresponding to the solution ρ = D/d, S = ρ in − ρ, R = L = P = 0, which is stable for (iii) the SP phase which is stable for γ C < βk L < γ C where and (iv) the RSP phase, whose stability we determine numerically. Since the dilution rate controls the overall growth rate and therefore the relative fitness differences between phenotypes, changing the dilution rate will impact the stability of the fixed points and can lead to removal or enhancement of regions that admit periodic dynamics in the phase diagram. Supplementary Figure 11 shows the representative phase diagrams of stable fixed points in the chemostat model for preventative and immune defenses, and in the low binding affinity or Lotka-Volterra limit where the infection rates are k I = k LV P with k LV = αn r /K m .

Supplementary Note 2. PHASE DIAGRAMS WITHOUT RESISTANCE
We consider the phage-bacteria system given by Eq. (12) in the main text, and set f R = 0 to model the system where the host has not acquired resistance to phage. Supplementary Fig. 4 shows the phase diagrams displaying the regions of stable S, SP and E fixed points in the absence of resistant phenotype. The eigenvalues at the host which make extinction stable for α > d/n r . The S phase {f S , f I } = {1, 0} is stable for βk L < γ. Thus, there is a bistable region at α > d/n r and βk L < γ which is shown in Supplementary Fig. 4. The SP phase is stable for α < d/n r and βk L > γ. In this reduced model, in the region where E is stable there can exist a second extinction fixed point, such that f S = 0 and f I = const > 0, and in which infected host and phage decay exponentially. This pathological fixed point is never stable in the full model with resistance, but in the reduced model we consider it part of the E phases.

Supplementary Note 3. SWITCHING MAINTAINS MULTIPLE PHAGE TYPES IN THE PREVENTATIVE DEFENSE MODEL
We consider two phage types P 1 and P 2 that compete for the same host. In this model, a host cannot be simultaneously infected by more than one phage type, but the infected cell can absorb phage of any type. The two phage types can have different absorption rates α and burst sizes β. The equations arė where A = S + I 1 + I 2 . We consider the minimal sensitivity condition n r = 1, but the model is easily generalizable for any n r > 1. Supplementary Fig. 5 shows the phase diagrams for the two-phage system, which include the RSP 1 , RSP 2 , and RSP 1 P 2 phases, where both phages stably coexist. The phase diagram is shown for P 2 given P 1 with coordinates {α 1 , β 1 } at three different locations in the RSP phase. For zero switching rate there will be no phage coexistence, since the RSP phase collapses to pure R, which does not support phage growth. Dashed curves denote transitions from single to multiple phage phases where both phage types are maintained in the RSP phase. In the region where P 2 outcompetes P 1 resistance can be lost if P 2 parameters are located in P 2 's SP phase; this region occurs at high P 2 burst rates and is shown in pink. The size and shape of these regions does not strongly depend on the value of s > 0.

Supplementary Note 4. FIXED POINT ANALYSIS FOR PATCH INVASION DYNAMICS
We analyze the replicator dynamics for the payoff matrix φ given in the main text, which we rewrite as where a 1 ≡ 2g 1 /(g 1 + c), a 2 ≡ 2g 2 /(g 2 + c), and a 3 ≡ 2s/(s + c). We always have 0 ≤ a i < 2, and additionally we assume a 1 ,a 2 > 1. Supplementary Table I lists the possible fixed points, which include the three pure strategies and two mixed strategies, y 1 and y 2 . y 1 is a fixed point on the R s -R 0 edge for a 3 < 1 (i.e. s < c), and a Nash equilibrium if Sφy 1 ≤ y 1 φy 1 which occurs for a 3 ≥ 1 − 1/a 2 . It is not an ESS because y 1 does not beat alternative best replies (e.g. R 0 φy 1 = y 1 φy 1 , yet y 1 φR 0 = (1 − a 3 )/(2 − a 3 ) < 1 = R 0 φR 0 ). y 2 is a Nash equilibrium in the simplex interior, hence each of the pure strategies is an alternative best reply. For y 2 to be an ESS, it must beat each of the pure strategies, i.e. y 2 φx > xφx must hold for x = R s , R 0 , and S. This implies a 1 − 1 > r 2 /r 1 , a 3 − 1 > r 3 /r 2 , and a 2 − 1 > r 1 /r 3 . Multiplying the inequalities we obtain (a 1 − 1)(a 2 − 1)(a 3 − 1) > 1 which cannot hold since 0 ≤ a i < 2, hence y 2 is not an ESS.
For a 3 ≥ 1 (s ≥ c), there are no stable fixed points, the interior fixed point y 2 is an unstable focus, and flows in the interior approach the heteroclinic cycle on the boundary. For decreasing a 3 , R s becomes a Nash equilibrium at a 3 = 1 (s = c), and then bifurcates to yield the Nash equilibrium y 1 and the unique ESS R s . As a 3 decreases further, y 2 moves toward the R s -R 0 edge, and at a 3 = (a 2 − 1)/a 2 it merges with y 1 and the latter loses its Nash equilibrium status. At each bifurcation a qualitative change occurs in the flow structure (see Supplementary Fig. 12 for representative plots).

Supplementary Note 5. INDEPENDENT DISPERSION OF PHAGE
To model phage that can migrate independently of bacteria, we add the strategy P to the game, which corresponds to phage without host cells. Co-invasion events can then involve P and one of the three bacterial strains R s , R 0 , or S. The game R s vs. P resolves immediately with R s , as phage is already present with the R s strain. Likewise, R 0 vs. P resolves in favor of R 0 because phage cannot grow on resistant cells and therefore cannot establish. Lastly, in S vs. P, phage can grow on sensitive cells, and drives S to extinction. We assign a 1 as the payoff to P against S, and note that R s outcompetes S by the same mechanism, i.e. phage-mediated killing of sensitive cells. Since R s additionally contains exponentially growing resistant cells, we have a 1 ≥ a 1 . The patch invasion game dynamics thus have the payoff matrix If a 1 < a 1 , P is a dominated strategy and is therefore eliminated, reducing the 4-strategy game given with Supplementary Eq. (9) to the R s R 0 S face discussed above. For a 1 = a 1 , the fixed points of the full 4-strategy game are the same those listed in Table I with two additional fixed points: P and with Q ≡ (1 + a 1 )(1 + a 2 ). One can verify that the conditions given in Table I for Nash equilibria or ESS hold for the 4-strategy game, that neither P nor y 3 are Nash equilibria of Supplementary Eq. (9), and that there is no fixed point in the simplex interior.
For the reduced game on the R 0 SP face, y 3 is a Nash equilibrium, but it is not an ESS (by the same reasoning used above to analyze y 2 ). Trajectories in the interior of the R 0 SP face approach the heteroclinic cycle on the face boundary. Since y 3 is a Nash equilibrium only in the reduced game, the eigenvector of the Jacobian in the direction transverse to the R 0 SP face has a positive eigenvalue, i.e. it is unstable.
In summary, the long-term outcome for the 4-strategy game given with Supplementary Eq. (9) with independent dispersion of phage is similar to that of the 3-strategy game given with Supplementary Eq. (8) of the main text. For s > c, there exist no stable fixed points and we find that trajectories in the simplex interior approach the heteroclinic cycle on the R s R 0 S face, while for s < c, R s is the unique ESS of the system.