1. Computational and Systems Biology
  2. Evolutionary Biology
Download icon

Repeated outbreaks drive the evolution of bacteriophage communication

  1. Hilje M Doekes  Is a corresponding author
  2. Glenn A Mulder
  3. Rutger Hermsen
  1. Theoretical Biology, Department of Biology, Utrecht University, Netherlands
  2. Laboratory of Genetics, Department of Plant Sciences, Wageningen University, Netherlands
Research Article
  • Cited 1
  • Views 1,432
  • Annotations
Cite this article as: eLife 2021;10:e58410 doi: 10.7554/eLife.58410

Abstract

Recently, a small-molecule communication mechanism was discovered in a range of Bacillus-infecting bacteriophages, which these temperate phages use to inform their lysis-lysogeny decision. We present a mathematical model of the ecological and evolutionary dynamics of such viral communication and show that a communication strategy in which phages use the lytic cycle early in an outbreak (when susceptible host cells are abundant) but switch to the lysogenic cycle later (when susceptible cells become scarce) is favoured over a bet-hedging strategy in which cells are lysogenised with constant probability. However, such phage communication can evolve only if phage-bacteria populations are regularly perturbed away from their equilibrium state, so that acute outbreaks of phage infections in pools of susceptible cells continue to occur. Our model then predicts the selection of phages that switch infection strategy when half of the available susceptible cells have been infected.

eLife digest

Bacteriophages, or phages for short, are viruses that need to infect bacteria to multiply. Once inside a cell, phages follow one of two strategies. They either start to replicate quickly, killing the host in the process; or they lay dormant, their genetic material slowly duplicating as the bacterium divides. These two strategies are respectively known as a ‘lytic’ or a ‘lysogenic’ infection.

In 2017, scientists discovered that, during infection, some phages produce a signalling molecule that influences the strategy other phages will use. Generally, a high concentration of the signal triggers lysogenic infection, while a low level prompts the lytic type. However, it is still unclear what advantages this communication system brings to the viruses, and how it has evolved.

Here, Doekes et al. used a mathematical model to explore how communication changes as phages infect a population of bacteria, rigorously testing earlier theories. The simulations showed that early in an outbreak, when only a few cells have yet been infected, the signalling molecule levels are low: lytic infections are therefore triggered and the phages quickly multiply, killing their hosts in the process. This is an advantageous strategy since many bacteria are available for the viruses to prey on. Later on, as more phages are being produced and available bacteria become few and far between, the levels of the signalling molecule increase. The viruses then switch to lysogenic infections, which allows them to survive dormant, inside their host.

Doekes et al. also discovered that this communication system only evolves if phages regularly cause large outbreaks in new, uninfected bacterial populations. From there, the model was able to predict that phages switch from lytic to lysogenic infections when about half the available bacteria have been infected.

As antibiotic resistance rises around the globe, phages are increasingly considered as a new way to fight off harmful bacteria. Deciphering the way these viruses communicate could help to understand how they could be harnessed to control the spread of bacteria.

Introduction

For several decades now, it has been recognised that communication between individuals is not limited to multicellular organisms, but is also common among microbes. The best-known example of microbial communication is bacterial quorum sensing, a process in which bacteria secrete signalling molecules to infer the local cell density and consequently coordinate the expression of certain genes (Nealson et al., 1970; Miller and Bassler, 2001). A wide variety of bacterial behaviours have been found to be under quorum-sensing control (Miller and Bassler, 2001; Hense and Schuster, 2015), including bioluminescence (Nealson et al., 1970), virulence (Antunes et al., 2010), cooperative public good production (Diggle et al., 2007; Darch et al., 2012), and antimicrobial toxin production (Cornforth and Foster, 2013; Kleerebezem and Quadri, 2001). Remarkably, it has recently been discovered that even some bacterial viruses (bacteriophages or phages for short) use signalling molecules to communicate (Erez et al., 2017). Here, we use a mathematical model to explore the dynamics of this viral small-molecule communication system. We study under what conditions communication between phages evolves and predict which communication strategies are then selected.

Bacteriophages of the SPbeta group, a genus in the order of Caudovirales of viruses that infect Bacillus bacteria, encode a small signalling peptide, named ‘arbitrium’, which is secreted when the phages infect bacteria (Erez et al., 2017). These phages are temperate viruses, meaning that each time a phage infects a bacterium, it makes a life-cycle decision: to enter either (i) the lytic cycle, inducing an active infection in which tens to thousands of new phage particles are produced and released through host-cell lysis, or (ii) the lysogenic cycle, inducing a latent infection in which the phage DNA is integrated in the host cell’s genome (or episomally maintained) and the phage remains dormant until it is reactivated. This lysis-lysogeny decision is informed by the arbitrium produced in nearby previous infections: extracellular arbitrium is taken up by cells and inhibits the phage’s lysogeny-inhibition factors, thus increasing the propenstiy towards lysogeny of subsequent infections (Erez et al., 2017). Hence, peptide communication is used to promote lysogeny when many infections have occurred. Similar arbitrium-like systems have now been found in a range of different phages (Stokar-Avihail et al., 2019). Notably, these phages each use a slightly different signalling peptide and do not seem to respond to the signals of other phages (Erez et al., 2017; Stokar-Avihail et al., 2019).

The discovery of phage-encoded signalling peptides raises the question of how this viral communication system evolved. While the arbitrium system has not yet been studied theoretically, previous work has considered the evolution of lysogeny and of other phage-phage interactions. Early modelling work found that lysogeny can evolve as a survival mechanism for phages to overcome periods in which the density of susceptible cells is too low to sustain a lytic infection (Stewart and Levin, 1984; Maslov and Sneppen, 2015). In line with these model predictions, a combination of modelling and experimental work showed that selection pressures on phage virulence change over the course of an epidemic, favouring a virulent phage strain early on, when the density of susceptible cells is high, but a less virulent (i.e. lysogenic) phage strain later in the epidemic, when susceptible cells have become scarce (Berngruber et al., 2013; Gandon, 2016). Other modelling work has shown that if phages, lysogenised cells, and susceptible cells coexist for long periods of time, less and less virulent phages are selected (Mittler, 1996; Wahl et al., 2018). This happens because phage exploitation leads to a low susceptible cell density, and hence a virulent strategy in which phages rapidly lyse their host cell to release new phage particles that can then infect other cells no longer pays off (because few cells are available to infect).

Erez et al., 2017 propose that the arbitrium system may have evolved to allow phages to cope with the changing environment during an epidemic, allowing the phages to exploit available susceptible bacteria through the lytic cycle when few infections have so far taken place and hence the concentration of arbitrium is low, while entering the lysogenic cycle when many infections have taken place and the arbitrium concentration has hence increased. This explanation resembles results for other forms of phage-phage interaction previously found in Escherichia coli-infecting phages (Abedon, 2017; Abedon, 2019). In the obligately lytic T-even phages, both the length of the latent period of an infection and the subsequent burst size increase if additional phages adsorb to the cell while it is infected – a process called lysis inhibition (Hershey, 1946; Doermann, 1948; Abedon, 2019). In the temperate phage λ, the propensity towards lysogeny increases with the number of co-infecting virions, called the multiplicity of infection (MOI) (Kourilsky, 1973). In both cases, modelling work has shown that the effect of the number of phage adsorptions on an infection can be selected as a phage adaptation to host-cell density, as it allows phages to switch from a virulent infection strategy (i.e. a short latent period or a low lysogeny propensity) when the phage:host-cell ratio is low to a less virulent strategy (i.e. a longer latent period or higher lysogeny propensity) when the phage:host-cell ratio is high (Abedon, 1989; Abedon, 1990; Sinha et al., 2017).

Here, we present a mathematical model to test if similar arguments can explain the evolution of small-molecule communication between viruses, and to explore the ecological and evolutionary dynamics of temperate phage populations that use such communication systems. We show that arbitrium communication can indeed evolve and that communicating phages consistently outcompete phages with non-communicating bet-hedging strategies. We however find that communication evolves under certain conditions only, namely if the phages regularly cause new outbreaks in substantial pools of susceptible host cells. Moreover, when communication evolves under such conditions, we predict that a communication strategy is selected in which phages use arbitrium to switch from a fully lytic to a fully lysogenic strategy when approximately half of all susceptible cells have been infected. Finally, we investigate how reliable the arbitrium signal needs to be for such communication to evolve, and show that the results are remarkably robust against variation in the density of bacteria.

Materials and methods

Model

Request a detailed protocol

Following earlier models (e.g. Stewart and Levin, 1984; Berngruber et al., 2013; Sinha et al., 2017; Wahl et al., 2018), we use ordinary differential equations to describe a well-mixed system consisting of susceptible bacteria, lysogens (i.e. lysogenically infected bacteria), and free phages, but extend this system to include an arbitrium-like signalling peptide (Figure 1A). For simplicity, we consider phages that do not affect the growth of lysogenised host cells; susceptible bacteria and lysogens hence both grow logistically with the same growth rate r and carrying capacity K. Lysogens are spontaneously induced at a low rate α, after which they lyse and release a burst of B free phages per lysing cell. Free phage particles decay at a rate δ and adsorb to bacteria at a rate a. Adsorptions to lysogens result in the decay of the infecting phage, thus describing the well-known effect of superinfection immunity (Hutchison and Sinsheimer, 1971; Susskind et al., 1974; McAllister and Barrett, 1977; Kliem and Dreiseikelmann, 1989; Bondy-Denomy et al., 2016), whereas adsorptions to susceptible bacteria result in infections with success probability b. We consider the lytic cycle to be fast compared to both bacterial growth and the lysogenic cycle (Stewart and Levin, 1984; Berngruber et al., 2013; Sinha et al., 2017; Wahl et al., 2018), so that a lytic infection can be modelled as immediate lysis releasing a burst of B free phages. Since the genes encoding arbitrium production are among the first genes to be expressed when a phage infects a host cell (Erez et al., 2017; Stokar-Avihail et al., 2019), each infection leads to an immediate increase of the arbitrium concentration A by an increment c. The lysis-lysogeny decision is effected by the current arbitrium concentration: a fraction φ(A) of the infections results in the production of a lysogen, while the remaining fraction (1-φ(A)) results in a lytic infection. Arbitrium does not decay spontaneously in the model (since it is a small peptide, spontaneous extracellular degradation is considered to be negligible), but it is taken up by bacteria at a rate u (e.g. through general bacterial peptide importers such as OPP [Erez et al., 2017]), and then degraded intracellularly, thus reducing the arbitrium concentration A.

Model overview.

(A) Free phages infect susceptible bacteria, at which point a fixed amount of arbitrium is produced. This arbitrium is taken up and degraded by susceptible cells and lysogens. Upon infection, a cell enters the lysogenic cycle with propensity φ(A), or the lytic cycle with propensity (1-φ(A)); the lysogeny propensity φ(A) depends on the current arbitrium concentration. The lytic cycle leads to immediate lysis of the host cell and release of a burst of new virions. In the lysogenic cycle, the phage remains dormant in the lysogen population, which grows logistically with the same rate as the susceptible cell population. Lysogens are spontaneously induced at a low rate, at which point they re-enter the lytic cycle. (B) In communicating phages, the lysogeny propensity φ(A) is modelled by a step-function characterised by two phage characteristics: θ, the arbitrium concentration above which the phage increases its lysogeny propensity, and ϕmax, the lysogeny propensity of the phage at high arbitrium concentration.

Consider competing phage variants i that differ in their (arbitrium-dependent) lysogeny propensity φi(A). The population densities (cells or phages per volume unit) of susceptible bacteria S, phage particles Pi and corresponding lysogens Li, and the concentration of arbitrium A can then be described by:

(1)dSdt=rS(1N/K)logistic growthbaSiPiinfection,(2)dLidt=rLi(1N/K)logistic growth+φi(A)baSPilysogenic infectionαLiinduction,(3)dPidt=BαLiburst from induction+B(1φi(A))baSPiburst from lytic infectionδPiphage decayaNPiadsorption,(4)dAdt=cbaSiPi production upon infectionuNAadsorption and degradation,

where N=S+iLi is the total density of bacteria.

We study two scenarios for the lysis-lysogeny decision: (i) a baseline scenario in which the arbitrium concentration does not affect the lysis-lysogeny decision; each phage variant has a constant lysogeny propensity ϕi and (ii) a full scenario in which the arbitrium concentration does affect the lysis-lysogeny decision; each phage variant causes lytic infection when the arbitrium concentration is low, but switches to some lysogeny propensity ϕmaxi when the arbitrium concentration exceeds the phage’s response threshold θi (Figure 1B; note that we use ϕmaxi to denote a constant characteristic of the phage and φi to denote the function describing how phage variant i’s lysogeny propensity depends on the arbitrium concentration). In this second scenario, phage variants with constant lysogeny propensity are still included: variants with a response threshold θi=0 cause lysogenic infections with lysogeny propensity ϕmaxi independent of the arbitrium concentration. Scenario (ii) is hence an extension of scenario (i).

On top of the ecological processes described in Equation 1–4, the model also allows for evolution of the phages due to mutations that change the characteristics ϕ(max) and θ of phage variants. In Equation 1–4 terms describing these mutations were omitted for readability; they are described in detail in Appendix A1.1. In short, replication of any phage variant was assumed to produce mutants with slightly different characteristics (e.g. a slightly higher or lower lysogeny propensity) with a small probability μ. Under scenario (ii), mutations changing ϕmax and θ are implemented as independent processes.

Serial passaging

Request a detailed protocol

In natural settings as well as in some laboratory experiments, phages regularly cause large outbreaks in pools of susceptible cells that were previously unavailable to the phages (e.g. when phages are spread to a new area, or when phages are serially passaged in a lab setting). Such outbreaks perturb the phage and cell populations away from their equilibrium. To mimic such repeated perturbations, we expose the system of Equation 1–4 to a phage serial-passaging regime (mimicking the experimental set-up of, for example, Bull et al., 1993; Bull et al., 2004; Bollback and Huelsenbeck, 2007; Betts et al., 2013; Broniewski et al., 2020). We initialise the model with a susceptible bacterial population at carrying capacity (S=K cells per mL) and a small phage population (iPi=106 phages per mL) and numerically integrate Equation 1–4 for a time of T hours. Then a fraction of the phage population is taken and transferred to a new population of susceptible bacteria at carrying capacity and Equation 1–4 are again integrated for T hours. This cycle is repeated to bring about a long series of epidemics. Throughout the manuscript, a dilution factor of D=0.01 is used (i.e. the passaged sample is 1% of the phage population). Passaging does not alter the relative frequency of the different phage variants, thus ensuring that the phage variants that were highly prevalent in the phage population at the end of an episode remain at a high relative frequency at the start of the new episode.

In this set-up, only phages are passaged from one epidemic episode to the next. To assess the robustness of simulation results to changes in this protocol, a second set-up was considered in which a fraction of the full sample (susceptible cells, lysogens, phages, and arbitrium) was passaged. We furthermore tested how the results are affected by variation in the bacterial carrying capacity. For this, at the start of each episode a carrying capacity value was sampled from a gamma distribution with mean K. We control the level of noise through the variance of this gamma distribution.

Parameters

In total, the model (Equation 1–4) has nine parameters (excluding the phage characteristics ϕi, ϕmaxi, and θi, which vary between phage variants present in any given simulation). As far as we are aware, none of these have been estimated for phages of the SPBeta group, but many have been measured for other phages, most of which infect E. coli (Table 1, estimates taken from Little et al., 1999; De Paepe and Taddei, 2006; Wang, 2006; Shao and Wang, 2008; Zong et al., 2010; Berngruber et al., 2013). To reduce the number of parameters in our analysis, we nondimensionalised the equations to obtain five scaled parameter values (Appendix A1.3) and used the literature estimates to derive default values for these scaled parameters (Table 1). To account for the uncertainty in these estimates, we performed parameter sweeps consisting of 500 simulations with parameter values randomly sampled from broad parameter ranges (Table 1). To ensure that low values of the parameters were well-represented, parameter values were sampled log-uniformly.

Table 1
Model parameters.
Original parameters
ParameterDescription (dimension)Literature estimatesReferences
rNet replication rate of bacteria (hour−1)1.0Berngruber et al., 2013
KCarrying capacity of bacteria (cells mL−1)109Berngruber et al., 2013
aAdsorption rate of phages to bacteria (hour−1 (cells per mL)−1)10-910-7De Paepe and Taddei, 2006; Shao and Wang, 2008
bProportion of adsorptions of a phage to a susceptible cell that leads to infection (cells phage−1)set at 10-2, not measuredBerngruber et al., 2013
BBurst size (phages)10—3.5· 103De Paepe and Taddei, 2006; Wang, 2006
αRate of spontaneous lysogen induction (hour−1)10-4—10-3Little et al., 1999; Zong et al., 2010; Berngruber et al., 2013
δSpontaneous decay rate of free phages (hour−1)10-3—2·10-2De Paepe and Taddei, 2006
uUptake rate of arbitrium by cells (arbitrium mL−1 (cells per mL)−1)no estimates known-
Scaled dimensionless parameters used in parameter sweeps
ParameterDescriptionDefault valueParameter sweep range
B^=bBEffective burst size21—103
a^=aKrScaled adsorption rate of phages to cells101—100
δ^=δrScaled decay rate of phage particles0.0110-3—0.1
α^=αrScaled spontaneous phage induction rate10-310-4—10-2
u^=uKrScaled rate of arbitrium uptake and degradation by cells0.110-3—1
DDilution factor of phages at serial passages0.0110-3—0.1

Model analysis

Request a detailed protocol

Numerical integration was performed in Matlab R2017b, using the default built-in ODE-solver ode45. Scripts are available from https://github.com/hiljedoekes/PhageCom.

Next to numerical integration results, we analytically found expressions for the model equilibria and derived expressions for the evolutionarily stable strategy (ESS) under serial passaging in both scenarios (excluding and including arbitrium communication). Detailed derivations are provided in Appendix A3.

Results

Evolution of the lysis-lysogeny decision and arbitrium communication requires perturbations away from equilibrium

A common approach to analysing ODE-models such as Equation 1–4 is to characterise the model’s equilibrium states (Stewart and Levin, 1984; Wahl et al., 2018; Cortes et al., 2019). Such an analysis is provided in Appendix A2. However, we will here argue that to understand the evolution of arbitrium communication, and the lysis-lysogeny decision in general, considering the equilibrium states is insufficient.

Firstly, the function of the arbitrium system is to allow phages to respond to changes in the density of susceptible cells and phages as reflected in the arbitrium concentration. But when the system approaches an equilibrium state, the densities of susceptible cells and phages become constant, and so does the arbitrium concentration. Equilibrium conditions hence defeat the purpose of small-molecule communication such as the arbitrium system. Evolution of small-molecule communication must be driven by dynamical ecological processes, and hence can only be studied in populations that are regularly perturbed away from their ecological steady state.

Secondly, under equilibrium conditions natural selection can act on the lysis-lysogeny decision only if infections still take place, and hence lysis-lysogeny decisions are still taken. We argue that this is unlikely. If the phage population is viable (i.e. if the parameter values are such that the phages proliferate when introduced into a fully susceptible host population), the model converges to one of two qualitatively different equilibria, depending on parameter conditions (Appendix A2): either (i) susceptible host cells, lysogens and free phages all coexist, or (ii) all susceptible host cells have been infected so that only lysogens and free phages remain. The evolution of a constant lysogeny propensity in a host-phage population with a stable equilibrium of type (i) was recently addressed by Wahl et al., who show that under these conditions selection always favours phage variants with high lysogeny propensity (i.e. ϕ=1) (Wahl et al., 2018). However, only a narrow sliver of parameter conditions permits a stable equilibrium of type (i) (Sinha et al., 2017; Cortes et al., 2019), and when we estimated reasonable parameter conditions based on a variety of well-studied phages, we found that they typically lead to a stable equilibrium of type (ii) (Appendix A2, parameter estimates based on Little et al., 1999; De Paepe and Taddei, 2006; Wang, 2006; Shao and Wang, 2008; Zong et al., 2010; Berngruber et al., 2013). This is because phage infections tend to be highly effective: their large burst size and consequent high infectivity cause temperate phages to completely deplete susceptible host cell populations, replacing them with lysogens that are immune to superinfection and hence have a strong competitive advantage over the susceptible cells (Bossi et al., 2003; Gama et al., 2013). Under common parameter conditions, after a short epidemic the susceptible cells population is hence depleted and no more infections take place, causing the competition between different phage variants to cease (see Figure 2A for example dynamics). Then there is no long-term selection on the lysis-lysogeny decision, and studying its evolution in this state is pointless.

Figure 2 with 2 supplements see all
Results in the absence of phage communication.

(A) Short-term model dynamics under default parameter conditions (Table 1) and a passaging episode duration of T=12 h. The model was initialised with a susceptible bacterial population at carrying capacity (S=K) and a low frequency of phages (iPi=10-5KB), and upon passaging the phages were diluted 100-fold. Phage variants differ in their lysogeny propensity ϕi. Dynamics within a single passaging episode are further illustrated in Figure 2—figure supplement 1. (B) Long-term model dynamics for default parameter settings and T=12 h. Over many passages, a single phage variant (ϕ=0.04) is selected. (C) Distribution of phage variants at evolutionarily steady state as a function of the time between passages, T. A total number of 101 phage variants was included, with lysogeny propensities varying between ϕ1=0 and ϕ101=0.5. When the interval between passages is short, the susceptible cells are not depleted during the rounds of infection and a fully lytic strategy (ϕ=0) is selected. For larger values of T, however, a bet-hedging strategy with small but non-zero ϕ-value is selected (ϕ0.04). Almost identical results are obtained if the serial passaging set-up is altered to simulate serial passaging of a full sample (phages, susceptible bacteria, and lysogens) instead of phages only (Figure 2—figure supplement 2). (D) Parameter sweep results. The model was run with 500 sets of randomly sampled parameters, and for each run the most abundant ϕ-value in the population at evolutionarily steady state was plotted against the analytically predicted evolutionarily stable strategy (ESS; see Appendix A3 and Box 1). The dotted line is the identity line. The analytically derived ESS is a good predictor of the simulation outcome.

We therefore consider a scenario in which the phage and cell populations are regularly perturbed away from equilibrium. To do so, we simulate serial-passaging experiments by periodically transferring a small fraction of the phages to a new population of susceptible host cells at carrying capacity, thus simulating cycles of repeated outbreaks (see Materials and methods).

In the absence of arbitrium communication, bet-hedging phages are selected with low constant lysogeny propensity

To form a baseline expectation of the evolution of the lysis-lysogeny decision under the serial-passaging regime, we first considered a population of phage variants that do not engage in arbitrium communication, but do differ in their constant lysogeny propensity ϕi. Under typical parameter conditions (default values in Table 1), each passaging episode starts with an epidemic in which the susceptible cell population is depleted, followed by a period in which the bacterial population is made up of lysogens only (Figure 2A, dynamics shown for a passaging episode length T=12 h). The composition of the phage and lysogen populations initially changes over subsequent passaging episodes (Figure 2A), but eventually an evolutionarily steady state is reached in which one phage variant dominates the phage population (ϕ=0.04; Figure 2B), confirming that the lysis-lysogeny decision is indeed under selection.

The distribution of phage variants at evolutionarily steady state depends on the time between passages, T (Figure 2C). If this time is short (T5 h), the phage variant with ϕ=0 dominates at evolutionarily steady state. This is an intuitive result: under these conditions phages are mostly exposed to environments with a high density of susceptible cells, in which a lytic strategy is favourable. Surprisingly, however, if the time between passages is sufficiently long (T>5 h), the viral population at evolutionarily steady state always centres around the same phage variant, independent of T (ϕ=0.04; Figure 2C). This result can be explained by considering the dynamics within a passaging episode (see Figure 2A): Once the susceptible cell population has collapsed, free phages no longer cause new infections and are hence ‘dead ends’. New phage particles are then formed by reactivation of lysogens only, so that the distribution of variants among the free phages comes to reflect the relative variant frequencies in the lysogen population. Hence, when the time between passages is sufficiently long, the phage type that is most frequent in the sample that is eventually passaged is the one that is most frequent in the population of lysogens (Figure 2—figure supplement 1). Under default parameter conditions, these are the phages with a low lysogeny propensity of ϕ=0.04.

Note that although a single phage variant clearly dominates the population, some diversity is maintained (Figure 2B–C). This is due to a mutation-selection balance: because mutants with slightly different ϕ-values continuously arise from the dominant phage variant and selection against these mutants is weak (due to their similarity to the dominant phage variant), the balance between influx of mutant variants by mutation and their efflux by selection results in the long-term presence of these mutants in the population. Such a population consisting of a dominant variant and its close mutants is called a quasi-species (Eigen, 1971).

Next, we assessed the robustness of the results to changes in the serial passaging protocol. In the standard protocol, only phages are passaged between episodes. If instead the passaged sample consists of the full system (susceptible cells, lysogens, and phages), almost identical results are obtained (Figure 2—figure supplement 2). This is again explained by realising that, as long as the time between passages is sufficiently long, the distribution of variants in the free phages is equal to the distribution in the lysogens. Since lysogens need to be induced to contribute to a new outbreak and the induction rate α is the same for all phage variants, the contribution of passaged lysogens to the new outbreak does not alter the relative frequency of phage variants.

To examine how these results depend on the model parameters, we determined which phage variant was most abundant at evolutionarily steady state for 500 randomly chosen parameter sets (see Table 1 for parameter ranges), always using a long time between passages (T=24 h). The selected ϕ-values for all parameter settings lie between ϕ=0 and ϕ=0.12 (y-axis of Figure 2D). We can hence conclude that selection favours phages with low but usually non-zero lysogeny propensities. These phages employ a bet-hedging strategy: throughout the epidemic they ‘invest’ a small part of their infection events in the production of lysogens, such that they are maximally represented in the eventual lysogen population.

To better understand how the lysogeny propensity ϕ that is selected depends on parameter values, we derived an analytical approximation for the evolutionarily stable strategy (ESS) under the serial-passaging regime if the time between passages is sufficiently long (Appendix A3.1–2). Because the phage dynamics during an epidemic affect the dynamics of the susceptible cells and vice versa, phage fitness is frequency dependent and the ESS is not found by a simple optimisation procedure, but by identifying the particular ϕ-value, denoted ϕ*, that maximises phage fitness given that this strategy ϕ* itself shapes the dynamics of the epidemic (Box 1). We find that the ESS can be approximated by the surprisingly simple expression

(5) ϕ*=1-(bB)-1log(BKP0),

where P0 is the density of phages at the start of a passaging episode. This approximation corresponds well with the results of the parameter sweep (Figure 2D), indicating that it indeed captures the most important factors shaping the evolution of the lysogeny propensity ϕ.

Box 1.

Lysogeny propensity of the evolutionarily stable strategy (ESS).

An evolutionarily stable strategy (ESS) is a strategy that cannot be invaded by any other strategy. In the context of the lysogeny propensity ϕ, it is the value ϕ* such that a population currently dominated by a phage with ϕ=ϕ* cannot be invaded by any phage variant with a different ϕ-value. A phage variant with ϕ=ϕi invading in a resident population with the same ϕ=ϕi always grows exactly like the resident. If this is the best possible invader, any other phage variant must perform worse than the resident and cannot invade. Hence, the ESS is the optimal response to itself. However, we still have to define what it means to be the ‘‘best possible invader’’ under the serial-passaging regime. Note that if the time between passages is sufficiently long, phages are selected on their ability to produce lysogens during the active epidemic (see Main Text). The optimal invader is hence the phage variant that, when introduced at a very low frequency, produces the most lysogens per capita between time t=0 and the time that the susceptible cell population collapses, TE. The ϕ-value of the optimal invader depends on TE (red line in plot): if the epidemic phase is short, lysogens have to be produced quickly and a high ϕ-value is optimal, while if the epidemic lasts longer, phages can profit more from lytic replication and a lower ϕ-value is optimal. In turn, however, the duration of the epidemic TE depends on the lysogeny propensity ϕ of the resident phage population (blue line in plot): phages with a lower value of ϕ replicate more rapidly and hence cause an earlier collapse of the susceptible population. The ESS is the value ϕ* that is optimal given the collapse time TE(ϕ*) that results when ϕ* itself is the resident strategy. Graphically, this value can be identified as the intersection of TE(ϕ) and ϕopt(TE) (the red and blue lines).

Box 1—figure 1
The ESS is found as the intersection of the curves of (i) the duration of the epidemic as a function of the lysogeny propensity of the resident (blue line), and (ii) the optimal lysogeny propensity of the invader given a fixed duration of the epidemic (red line).

Equation 5 shows that the ESS depends on the initial phage density in a passaging episode, P0, relative to the burst size B and maximal host-cell density K, and the effective burst size bB, which represents the expected number of progeny phages per phage that adsorbs to a susceptible bacterium. The ESS ϕ* decreases with the dilution factor of the phages upon passage (i.e. with lower P0). On the other hand, ϕ* increases with the effective burst size bB (note that (bB)-1 decreases when (bB) increases). Both effects can be intuitively understood by considering how these factors affect the duration of the epidemic, TE. If the phage density is low at the start of a passaging episode or if the phages have a small effective burst size, it takes a while before the phage population has grown sufficiently to cause the susceptible population to collapse. Since a lytic strategy is favoured early in the epidemic, when the susceptible cell density is still high, a longer epidemic favours phages with lower values of ϕ (see the red line in the figure in Box 1). On the other hand, if the initial phage density is high or if the phages have a high effective burst size, the susceptible cell population collapses quickly, phages have a much shorter window of opportunity for lysogen production and hence phages with higher ϕ-values are favoured.

If arbitrium communication is included, communicating phages are selected that switch from a fully lytic to a fully lysogenic strategy

Next, we included the possibility of arbitrium communication and let phage variants be characterised by two properties: their arbitrium response threshold, θi, and their lysogeny propensity when the arbitrium concentration exceeds their response threshold, ϕmaxi (see Figure 1B). We then again considered the dynamics of our model under a serial-passaging regime.

In Figure 3A, example dynamics are shown for three competing phage variants, all with ϕmax=1 but with different response thresholds θi. The arbitrium concentration increases over the course of the epidemic, and the phage variants switch from lytic infection to lysogen production at different times because of their different response thresholds.

Figure 3 with 2 supplements see all
Model dynamics if phage communication is included in the model.

(A) Short-term dynamics for default parameter conditions (Table 1) and the same serial-passaging regime as in Figure 2. This example shows the competition between three phage variants, all with ϕmax=1 but with varying response thresholds θ. (B) Distribution of phage variants at evolutionary steady state for varying passaging episode durations T. In total 441 phage variants were included in this analysis, covering all combinations of ϕmax between 0 and 1 and θ between 0 and cK with step sizes 0.05 and 0.05cK, respectively. When the interval between passages is very short, again a fully lytic strategy (ϕmax=0) is selected. For longer times between passages, however, we consistently see that a strategy with ϕmax=1 and θ0.65cK dominates the population. The results shown for T=2 h are representative for values of T4 h, while the results shown for T=12 h represent results obtained for T5 h (see Figure 3—figure supplement 1 for distributions for a large range of T-values). Almost identical results are obtained if instead of only phages a full sample (susceptible cells, lysogens, phages, and arbitrium) is passaged (Figure 3—figure supplement 2). (C) Rapid invasion by 'optimally' communicating phages into a population of phages with the 'optimal' bet-hedging strategy. The bet-hedging phages have ϕ=0.04 (see Figure 2C), while the communicating phages are characterised by ϕmax=1 and θ=0.66cK (see panel C). The communicating phage is initialised at a frequency of 1% of the bet-hedging phage.

Note that the maximum arbitrium concentration obtained during a passaging episode is approximately A=cK (Figure 3A). This is because during the epidemic the dynamics of the susceptible cell density are mostly determined by infection events and not so much by the (slower) bacterial growth. Since the arbitrium concentration increases by an increment c every time a susceptible cell is infected, the infection of all initial susceptible cells will result in an arbitrium concentration of A=cK (assuming that the degradation of arbitrium is also slow and can be ignored during the growth phase of the epidemic). The arbitrium concentration during the early epidemic then is a direct reflection of the fraction of susceptible cells that have so far been infected.

To study the evolution of arbitrium communication, we again considered the distribution of phage variants at evolutionary steady state for varying values of the time between passages, T. Similar to the results shown in Figure 2C, we find two main regimes (Figure 3B): if the time between passages is short (T<4 h, illustrated by T=2 h in Figure 3B), selection favours phage variants that only cause lytic infections (ϕmax=0); if the time between passages is sufficiently long (T5 h, illustrated by T=12 h in Figure 3B), the phage population is dominated by variants with ϕmax=1 and θ0.65cK. For 4T<5 h, we see a transition between these two regimes (Figure 3—figure supplement 1). If the time between passages is sufficiently long (T>5 h), phage variants are hence selected that switch from a completely lytic to a completely lysogenic strategy when the arbitrium concentration exceeds a certain threshold.

In the simulations of Figure 3B, phage variants could have emerged that use the bet-hedging strategy found in the absence of communication (in phage variants with θ=0, the lysogeny propensity is always ϕmax, independently of the arbitrium concentration), but this did not happen. We can hence conclude that any bet-hedging phage variants were outcompeted by variants that do use arbitrium communication. To underscore this conclusion, we simulated a competition experiment between the bet-hedging phage variant that was selected in the absence of communication and the communicating variant selected when arbitrium dynamics were included (Figure 3C). The communicating phage quickly invades on a population of bet-hedging phages and takes over, confirming that communication is indeed favoured over bet-hedging.

If a full sample (susceptible cells, lysogens, phages, and arbitrium) is passaged instead of phages only, again almost identical results are found (Figure 3—figure supplement 2). As was the case for the simulations in which arbitrium was absent, passaged lysogens do not alter the distribution of phage variants in the new outbreak. The passaged arbitrium does not significantly affect the outbreak dynamics either, because its concentration after dilution is much lower than the response threshold θ of the phage variants that are selected.

Evolved phages switch from the lytic to the lysogenic life-cycle when approximately half of the susceptible cells have been infected

To study how the evolution of phage communication depends on phage and bacterial characteristics, 500 simulations were performed with randomly sampled sets of parameter values (Table 1), using a long time between serial passages (T=24 h). For each simulation, we determined which phage variant was most prevalent at evolutionary steady state. Although we varied the parameter values over several orders of magnitude, the most prevalent phage variant had a lysogeny propensity of ϕmax=1 and a response threshold of θ=0.5cK or θ=0.6cK in almost all simulations (Figure 4A). Hence, over a broad range of parameter values, phages are selected that use the arbitrium system to switch from a fully lytic to a fully lysogenic strategy (i.e. ϕmax=1). This suggests that over the course of an epidemic, there is an initial phase during which the lytic strategy is a ‘better’ choice (i.e. produces the most progeny on the long run), while later in the epidemic the production of lysogens is favoured and residual lytic infections that would results from a lysogeny propensity φ<1 are selected against.

Parameter dependance of the selected values of ϕ𝐦𝐚𝐱 and θ.

(A) Parameter sweep results. A total of 500 simulations were run with randomly sampled parameters (Table 1) and a long time between passages (T=24 h). In each simulation, 121 phage variants were included, covering all combinations of ϕmax=0 to ϕmax=1 and θ=0 to θ=cK with step sizes 0.1 and 0.1cK, respectively. The size of the circles corresponds to the number of simulations that yielded that particular phage variant as most abundant at evolutionary steady state. (B) Analytically predicted θ-value as a function of the effective burst size per adsorption to a susceptible cell, bB, and most abundant phage variant found in a simulation with varying bB but otherwise default parameter values, T=24 h, ϕmax=1 and θ=0,0.02cK,,cK. The range on the x-axis is equal to the range sampled in the parameter sweep. The analytically derived evolutionarily stable θ* is a good prediction for the response threshold selected in the simulations, especially for phages with high effective burst size.

To better understand the intriguing consistency in θ-values found in the parameter sweep, we used a similar approach as before to analytically derive an approximation for the response threshold θ* of the evolutionarily stable strategy under the condition that the time between passages is long (Appendix A3.3). Again, we find a surprisingly simple expression for the ESS:

(6) θ*=cK2-(bB)-1.

Note that the expression in Equation 6 again depends on the effective burst size bB, which is an indicator of the phage’s infectivity. The evolutionarily stable response threshold θ* declines as the effective burst size increases, converging to a value of θ*=12cK for highly infective phages (Figure 4B, green line). The same result was found for simulations of the competition between phage variants with different θ-values under different effective burst sizes (Figure 4B, blue dots). We see that Equation 6 provides a good prediction for the response threshold value that is selected over evolutionary time, especially for phages with high effective burst size (Figure 4B).

For phages with a very small effective burst size, the response threshold selected in the simulations tends to be lower than the analytical approximation. This is due to a violation of one of the simplifying assumptions made to arrive at the analytical approximation of Equation 6, namely that during the active epidemic the dynamics of the arbitrium concentration are dominated by its production through infections and arbitrium uptake and degradation by susceptible cells can be ignored. While this is a reasonable assumption in a fast progressing epidemic, it breaks down if the dynamics of the epidemic are slow, which is exactly the case if the effective burst size bB is small. Under these conditions, the uptake and degradation of arbitrium by susceptible cells cause the arbitrium concentration to be lower than assumed in the analytical derivation. Consequently, the actual selected response thresholds (which are essentially arbitrium concentration values) are lower than the analytically predicted values.

The result in Equation 6 can be further understood biologically. Remember that the arbitrium concentration during the epidemic varies between A=0 and A=cK, and is a reflection of the fraction of susceptible cells that have so far been infected. It makes sense that the evolutionarily stable response threshold causes phages to switch infection strategy somewhere in the middle of the epidemic: if a phage variant switches to the lysogenic strategy too early, its free phage population does not expand enough to compete with phages that switch later, but if it switches too late, the susceptible-cell density has decreased to such a degree that the phage has missed the window of opportunity for lysogen production. The ESS results from a balance between the fast production of phage progeny during the initial lytic cycles and the eventual production of sufficient lysogens. For phages with a high effective burst size, this balance occurs around the time that half of the available susceptible cells have been infected. Phages with lower effective burst size are, however, predicted to switch later, because these phages need to invest a larger portion of the available susceptible cells in the production of free phages to produce a sufficient pool of phages that can later form lysogens. Note, however, that the range of biologically reasonable effective burst sizes includes many high values (range of x-axis in Figure 4B, Table 1), that is, many real-life phages have high infectivity. Hence, for natural phages in general, we predict that if they evolve an arbitrium-like communication system, communication will be used to switch from causing mostly lytic to mostly lysogenic infections when in an outbreak approximately half of the pool of susceptible bacteria has been infected.

Arbitrium communication is robust against variation in bacterial carrying capacity

So far, we have considered the evolution of arbitrium communication under highly predictable settings, with each outbreak taking place in a population of bacteria with the same initial density (i.e. the bacterial carrying capacity was constant). As argued above, in such a set-up the arbitrium concentration provides information on the density of susceptible cells still available for infection, and the phages use this to inform their lysis-lysogeny decision. While the bacterial carrying capacity can be kept constant in lab experiments, it is far from obvious that this would be the case in natural environments. This warrants the question of how robust the results are to variation in the bacterial carrying capacity.

We therefore performed simulations in which the carrying capacity varies from outbreak to outbreak. For a long time between passages (T = 24 h), at the start of each passaging episode a random carrying capacity was drawn from a gamma distribution with mean K and a pre-set variance that differs from simulation to simulation. We use the coefficient of variation (CV), which is defined as the standard deviation relative to the mean, to describe the level of noise.

Figure 5 summarises the results of these additional simulations. Surprisingly, a communication strategy with ϕmax=1 and θ0.5cK is selected for a large range of carrying capacity noise up to CV ≤0.35 (illustrated by CV = 0.22 in Figure 5A; see Figure 5—figure supplement 1 for full data). In other words, even if the carrying capacity varies with a standard deviation up to one third of its mean value, the communication strategy described in the previous section is still selected.

Figure 5 with 1 supplement see all
Distribution of phage variants at evolutionary steady state for increasing variation in the bacterial carrying capacity.

Each simulation included a total of 441 phage variants that covered all combinations of ϕmax between 0 and 1 and θ between 0 and cK with step sizes 0.05 and 0.05cK, respectively. In each simulation, 1000 passaging episodes were simulated, with a long time between passages (T=24 h). Parameter were set to default values (Table 1), except that at the start of each passaging episode the value of the bacterial carrying capacity was drawn from a gamma distribution with mean K. The coefficient of variation (CV = standard deviation / mean) was varied between simulations. Results are shown here for (A) CV = 0.22, (B) CV = 0.39, (C) CV = 0.59. The results in panel (A) are representative for CV ≤ 0.35, and the results in panel (C) are representative for CV ≥ 0.5 (Figure 5—figure supplement 1).

As the coefficient of variation increases even further, the arbitrium response threshold value θ of the selected phages decreases, and so does the lysis-lysogeny propensity that is used at high arbitrium concentration ϕmax (Figure 5B and C). These results make sense: if the carrying capacity strongly varies between passaging episodes, the phages regularly cause outbreaks in bacterial populations with low density. Phages with a response threshold value larger than the bacterial carrying capacity do not produce any lysogens during such an outbreak, which is disastrous for their long-term fitness. Hence, lower response thresholds are selected. The corresponding lower ϕmax values likely evolve to compensate for the earlier switch to lysogen production caused by the lower θ-values. In highly variable conditions, phages are hence selected to switch from a lytic strategy very early in the epidemic to a bet-hedging strategy later.

While we find much lower response threshold values when the variation in bacterial carrying capacity is high, these threshold values do remain clearly larger than zero (Figure 5C). This is true even if the carrying capacity is exponentially distributed (CV = 1; see Figure 5—figure supplement 1). Hence, even under very high variation of bacterial density a form of arbitrium communication (in which phages use the arbitrium signal to switch from a lytic to a bet-hedging strategy) is still favoured over completely bet-hedging strategies.

Discussion

We have presented a mathematical model of a population of phages that use an arbitrium-like communication system, and used this model to explore the evolution of the lysis-lysogeny decision and arbitrium communication under a serial-passaging regime. When arbitrium communication was excluded from the model, we found that bet-hedging phages with relatively low lysogeny propensity were selected. But when arbitrium communication was allowed to evolve these bet-hedging phages were outcompeted by communicating phages. These communicating phages switch from a lytic strategy early in the epidemic to a fully lysogenic strategy when approximately half of the available susceptible cells have been infected.

The serial-passaging set-up of the model is crucial for the evolution of the lysis-lysogeny decision and arbitrium communication. This has two main reasons. Firstly, it ensures that the phages are regularly exposed to susceptible cells, thus maintaining selection pressure on the lysis-lysogeny decision. Because of their high infectivity (see Materials and methods section and De Paepe and Taddei, 2006; Wang, 2006), most temperate phage outbreaks will completely deplete pools of susceptible bacteria, resulting in a bacterial population consisting of lysogens only in which the phage no longer replicates through infection (Bossi et al., 2003; Gama et al., 2013). The bet-hedging strategy we found in the absence of phage communication is a mechanism to deal with these (self-inflicted) periods of low susceptible cell availability, consistent with earlier studies (Maslov and Sneppen, 2015; Sinha et al., 2017). Secondly, the serial-passaging set-up imposes a dynamic of repeated epidemics in which a small number of phages is introduced into a relatively large pool of susceptible cells. Such dynamics are necessary for the arbitrium system to function: the arbitrium concentration provides a reliable cue for a phage’s lysis-lysogeny decision only if it is low at the beginning of an epidemic and subsequently builds up to reflect the fraction of cells that have so far been infected.

Based on these considerations, we can stipulate which environments promote the evolution of small-molecule communication such as the arbitrium system. One major factor that can ensure a regular exposure to susceptible cells (the first requirement) is spatial structure. If phages mostly infect bacteria that are physically close to them, a global susceptible population can be maintained even though susceptible bacteria may be depleted in local environments (Kerr et al., 2006). Indeed, spatial structure has been shown to greatly influence phage evolution, for instance by promoting the selection of less virulent strains that deplete their local host populations more slowly (Kerr et al., 2006; Heilmann et al., 2010; Berngruber et al., 2015). For small-molecule communication to evolve, however, the phages would also have to undergo repeated, possibly localised, outbreak dynamics (the second requirement). Such dynamics could occur in structured meta-populations of isolated bacterial populations, between which the phages spread at a limited rate. Alternatively, phages might encounter large pools of newly susceptible bacteria if they escape superinfection immunity through mutation (Zinder, 1958; Bailone and Devoret, 1978; Scott et al., 1978). Under this scenario, however, any remaining arbitrium signal from previous infection events no longer provides accurate information about the number of susceptible cells available, since cells that were lysogenically infected are once again susceptible to infection with the new phage variant. If the escape mutation occurs after the arbitrium produced during previous epidemics has been degraded, this problem does not occur and the newly produced arbitrium does function as a reliable signal of susceptible cell density for the new phage variant. If, however, the escape mutation occurs while the arbitrium concentration is still high from previous outbreaks, the new phage variant will cause lysogenic infections while in fact the lytic cycle should be favoured. There will then be selection pressure on the new phage variant to acquire additional mutations that change its signal specificity. This might in part explain the large diversity of phage signalling peptides observed (Erez et al., 2017; Stokar-Avihail et al., 2019).

The model presented in this paper allows us to put hypotheses about the arbitrium system to the test. For instance, it has been suggested that the arbitrium system would benefit from the production of arbitrium by lysogens, because phages thereby would be warned about the presence of neighbouring lysogens (which are immune to superinfection; Hynes and Moineau, 2017). Above we have argued, however, that under repeated epidemics, such as caused by serial passaging, selection on the lysis-lysogeny decision and arbitrium signalling is limited to the relatively short window of time in which all (locally) present susceptible cells become infected: afterwards no new infections occur and arbitrium therefore has no effect. During this short time window, the density of lysogens is still low, and any arbitrium produced by lysogens contributes little to the information already conveyed by arbitrium produced during infection events. Hence, our model predicts that, under repeated epidemics that completely deplete (local) pools of susceptible cells, the effects of arbitrium production by lysogens are likely very minimal. Arbitrium production by lysogens can be effective only if lysogens and susceptible cells coexist over sufficiently long periods of time, such that infection events occur in the presence of lysogens. In the model, we found that such coexistence is highly unlikely. Coexistence between lysogens and susceptible cells might, however, happen under circumstances that were not included in our model, for instance through a constant inflow of susceptible cells because of cell migration, or through the loss of superinfection immunity by lysogens.

Intriguingly, our model predicts that phages using small-molecule communication to coordinate their lysis-lysogeny decision would be selected to switch from a lytic to a lysogenic strategy once approximately half of the available susceptible bacteria have been lytically infected. This prediction warrants experimental testing. However, it also raises the question of how the phages would ‘know’ at what bacterial density the susceptible population has been halved. For the arbitrium signal to carry reliable information about the density of remaining susceptible cells, the initial concentration of susceptible bacteria has to be similar from outbreak to outbreak. Hence, one might expect the communication strategy to break down if the density of susceptible bacteria is variable. Surprisingly, this turned out not to be the case. We found that arbitrium-like communication could evolve even if the bacterial carrying capacity was highly variable. The characteristics of the communication system then depend on the level of noise. In highly variable environments, we predict the selection of phages that start their lysogen production earlier in an outbreak (i.e. phages that have a low response threshold), and then do so in a bet-hedging way (i.e. with a lysogeny propensity much smaller than 1).

In fact, few details are known so far about the response curve of phages’ lysogeny propensity to the arbitrium concentration. In the model, we chose to implement the response to arbitrium as a stepwise function. This allowed us to clearly distinguish between strategies that are favoured at low arbitrium concentration (the lytic cycle) and at high arbitrium concentration (the lysogenic cycle). In reality, phages might respond more gradually to the arbitrium concentration. While this would alter some of our results (e.g. pinpointing an arbitrium concentration at which the phages switch infection strategy becomes harder, if not impossible), we do not expect the results in general to depend on the precise shape of the response curve: phages will still use the arbitrium signal to adjust their infection strategy to whichever strategy currently yields most progeny phage on the long run. Once more data become available on the actual shape of the response curve, these can be incorporated in the model by adjusting the arbitrium response function φ(A), thus producing a more specific model of the arbitrium system.

Next to the arbitrium system, several other examples of temperate phages affected by small signalling molecules have recently been described. For instance, the Vibrio cholerae-infecting phage VP882 ‘eavesdrops’ on a quorum-sensing signal produced by its host bacteria, favouring lytic over lysogenic infections when the host density is high (Silpe and Bassler, 2019), while in coliphages λ and T4 and several phages infecting Enterococcus faecalis, the induction of prophages, that is, the lysogeny-lysis decision, is affected by bacterial quorum sensing signals (Ghosh et al., 2009; Rossmann et al., 2015; Laganenka et al., 2019). The model could be adapted to capture these other regulation mechanisms by changing the arbitrium equation to an equation describing the production and degradation of the bacterial quorum sensing signal, and – for the second mechanism – letting the prophage reactivation rate α, rather than the lysogeny propensity ϕ, depend on the signal concentration. Similar analyses to the ones in this paper would then allow us to study under what conditions phage eavesdropping on bacterial quorum sensing can and cannot evolve. Mathematical and computational modelling can thus help to better understand the ecology and evolution of these fascinating regulation mechanisms as well.

Appendix 1

Model equations and parameters

A1.1 Full model equations, including mutations

For readability, the model equations in the main text (Equation 1–4) did not include mutations between phage variants. Here, we present the full model equations, including mutations. We assume that mutations happen when new phage particles (and hence new copies of the phage’s genetic material) are formed prior to lysis of the host cell. Infection with a parent phage of variant j results in progeny phages of variant i with probability μij, where μii is the probability that offspring of a parent phage of type-i has no mutations. The full model reads:

(A.1)dSdt=rS(1N/K)baSP,(A.2)dLidt=rLi(1N/K)+φi(A)baSPiαLi,(A.3)dPidt=Bjμij(αLj+[1φj(A)]baSPj)δPiaNPi,(A.4)dAdt=cbaSPuNA,(A.5)withφi(A)=ϕiif arbitrium communication is excluded from the model,(A.6)andφi(A)={0if Aθi,ϕmaxiif A>θiif arbitrium communication is included (see Figure 1B).

In the first part of the manuscript, where we exclude arbitrium communication from the model, phage variants are characterised by their constant lysogeny propensity ϕi (Equation A.5). Mutations between phage variants were implemented in a stepwise fashion, that is, μij=μ>0 if ϕj is one step higher or lower than ϕi, and μij=0(ij) otherwise. Throughout this study, a value of μ=510-4 was used. Varying μ alters the mutation-selection balance and thus affects the frequency of mutants in the quasi-species. As long as μ is reasonably small, however, it does not change which phage variant dominates the population, that is, μ does not affect the evolutionarily stable strategy.

In the second part of the manuscript, where we include the possibility of arbitrium communication, phage variants are characterised by two properties (Equation A.6): their arbitrium threshold θi and the lysogeny propensity obtained when the arbitrium concentration exceeds this threshold, ϕmaxi. Mutations in the values of ϕmax and θ were implemented as independent processes, both happening in a stepwise fashion as explained above.

A1.2 Phage variants included in simulations

In the simulations of the restricted model (arbitrium communication excluded) shown in Figure 2C, a range of phage variants was included with ϕ1=0, ϕ2=0.005, , ϕ100=0.495, ϕ101=0.5. In the simulations presented in Figure 2D, a range of phage variants was included of ϕ1=0, ϕ2=0.01, , ϕ100=0.99, ϕ101=1.0.

In the simulations of the full model (arbitrium communication included) shown in Figure 3B, 441 phage variants were included, representing all possible combinations of ϕmax=0,0.05,,1 and θ=0,0.05cK,,cK. In the simulations of Figure 4A, 121 phage variants were included representing all possible combinations of ϕmax=0,0.1,,1 and θ=0,0.1cK,,cK. In the simulations of Figure 4B, all phage variants had ϕmax=1, but they varied in θ=0,0.02cK,,cK.

All simulations were initialised with the susceptible cells at carrying capacity (S=K cells per mL), no lysogens (iLi=0) and a low total density of iPi=106 phages per mL. All phage variants were initially present at equal frequency.

In all simulations mutations between phage variants were included as described in the previous section, with the exception of the competition experiment between the optimal bet-hedging phage variant and the optimal communicating phage variant (results shown in Figure 3C). In this latter case, the mutation rate was set to zero.

A1.3 Parameter reduction

In total, the model of Equation A.1–A.6 has nine parameters (ignoring mutation probabilities). This number can, however, be reduced by non-dimensionalising the equations. We introduce the dimensionless variables

t^=rt,S^=SK,L^i=LiK,P^i=PiKB,A^=AcK,θ^i=θicK,

and define P^=iP^i, L^=iL^i, and N^=S^+L^. Let furthermore

B^=bB,a^=aKr,α^=αr,δ^=δr,u^=uKr.

Using these new variables and parameters, the equations reduce to:

(A.7)dS^dt^=S^(1N^)B^a^S^P^,(A.8)dL^idt^=L^i(1N^)+φ^i(A^)B^a^S^P^iα^L^i,(A.9)dP^idt^=jμij(α^L^j+[1φ^j(A^)]B^a^S^P^j)δ^P^ia^N^P^i,(A.10)dA^dt^=B^a^S^P^u^N^A^,(A.11)withφ^i(A^)=ϕiif arbitrium communication is excluded from the model,(A.12)andφi(A)={0if Aθ^i,ϕmaxiif A>θ^iif arbitrium communication is included.

Five dimensionless parameters are left in Equation A.7–A.12: the effective burst size per adsorption of a phage to a susceptible cell, B^, and the scaled rates a^, α^, δ^, and u^. These non-dimensionalised equations are used throughout the rest of this appendix, unless stated otherwise, dropping the hats for convenience.

Appendix 2

Equilibrium analysis

In this section, we find the dynamical equilibria existing in the model, and derive parameter conditions for their stability. This analysis provides us with a baseline expectation of the densities of phages, lysogens, and susceptible cells that the model converges to after sufficient time.

Equilibria of the model are found by equating Equation A.7–A.10 to zero and solving for the model variables. By the definition of equilibrium, the arbitrium concentration at equilibrium is constant:

(A.13) A¯={BaS¯P¯uN¯if N¯>0,A¯=0if N¯=0,

where the bar indicates equilibrium values. Because the equilibrium arbitrium concentration is constant, at equilibrium the different phage variants are characterised by their lysogeny propensity φi(A¯) only, irrespective of whether arbitrium communication is included in the model or not. Hence, expressions for the model equilibria are the same in the absence and presence of arbitrium communication. Below, we derive these expressions, and determine under what conditions the different model equilibria are stable.

In the model, four qualitatively different types of equilibria can occur.

Firstly, there is a trivial equilibrium at S¯=0, and L¯i=0, P¯i=0 for all phage variants i. This equilibrium is unstable as long as the bacterial logistic growth rate satisfies r>0.

Secondly, there is an equilibrium in which the susceptible population is at carrying capacity and the infection is absent: S¯=1, and L¯i=0, P¯i=0 for all phage variants i. This equilibrium is stable if no phage-lysogen pairs Pi-Li can invade on the susceptible population. To derive stability conditions, we approximate the dynamics of Pi and Li in the vicinity of the equilibrium by the linearised equations

(A.14) (dLidtdPidt)(-αϕ¯iBaα(1-ϕ¯i)Ba-δ-a)(LiPi)=:J(LiPi),

where ϕ¯i=φi(0) is the lysogeny propensity of phage type i at the equilibrium. No phage-lysogen pair can invade (i.e. the equilibrium is stable) precisely if the real parts of both eigenvalues of the Jacobian matrix J are negative for all i. The eigenvalues of the Jacobian matrix J are given by

λ+/=Γ±Δ2,with Γ=a(B(1ϕ¯i)1)δα and Δ=Γ2+4α(a(B1)δ).

The real parts of λ+/- are both negative precisely if Δ<Γ2 and Γ<0. From Γ<0, we find

(A.15) a(B(1-ϕ¯i)-1)<δ+α,

while Δ<Γ2 yields

(A.16) a(B-1)<δ.

Since all parameters are non-negative and 0ϕ¯i1, the condition in Equation A.16 is more stringent than the condition in Equation A.15. Note also that the condition in Equation A.16 does not depend on the lysogeny propensity of the invading phages, ϕ¯. Hence, the equilibrium S¯=1, L¯i=0, and P¯i=0 for all phage variants i is stable exactly if the condition in Equation A.16 is satisfied. This condition makes sense: phages cannot spread in a susceptible cell population at carrying capacity if their infection rate BaS¯ is smaller than the decay rate of phage particles δ+aS¯.

Thirdly, there is a class of equilibria in which S¯=0, and L¯i>0, P¯i>0 for some i. In these equilibria, the total densities of lysogens L¯ and of free phages P¯ are given by

(A.17) L¯=1-α,P¯=αL¯δ+aL¯=α(1-α)δ+a(1-α).

In the absence of susceptible cells at equilibrium, no infections can take place and hence all phage variants behave identically (since phages vary only in φi(A), which occurs exclusively in the infection terms). This is reflected in the equations for the different lysogen variants, which for S¯=0 and L¯=1-α (Equation A.17) reduce to

dLidt=Li(1-L¯)-αLi=0.

Hence, any combination of L¯i values with iL¯i=L¯=1-α and corresponding P¯i-values,

(A.18) P¯i=αjμijL¯jδ+a(1-α),

is an equilibrium. Analogous to the reasoning above, such an equilibrium is stable if the susceptible cells cannot invade the phage-lysogen population at equilibrium. The linearised equation for the dynamics of S near the equilibrium of Equation A.17 reads

(A.19) dSdtS(1-L¯)-BaSP¯=S(α-Baα(1-α)δ+a(1-α)).

The right-hand side of this equation is negative (i.e. susceptible cells cannot invade) precisely if δ+a(1-α)<Ba(1-α), or summarised

(A.20) δ<a(B-1)(1-α).

(Note that to arrive at condition Equation A.20 we assume δ+a(1-α)>0. This assumption is justified because the spontaneous induction rate of lysogens α is small, and hence 1-α>0 (see Table 1).)

Lastly, there can be an equilibrium in which the susceptible cells, lysogens and phages all coexist. Expressions for S¯, P¯i, and L¯i at this equilibrium are bulky and not directly insightful. This type of equilibrium was however extensively analysed in recent work by Wahl et al. for phages with a constant lysogeny propensity (i.e. the restricted model where arbitrium communication is excluded) (Wahl et al., 2018). Remember that in an equilibrium state, phage variants are characterised by their lysogeny propensity φi(A¯) only (i.e. differences in response threshold θ are relevant only if they are reflected in differences in φi(A¯); phage variants i and j with θiθj but φi(A¯)=φj(A¯) can for all practical purposes be considered the same), and hence the results found by Wahl et al. can be extended to the model analysed here. Phage variants with different lysogeny propensity φi(A¯) can be seen as consumers that compete for a single resource, namely susceptible cells to infect. As long as S¯>0, we hence expect competitive exclusion, and the phage population at equilibrium will be dominated by phages with a single lysogeny propensity value ϕ¯φ(A¯) (when mutations are ignored, these phages will be the only ones present; otherwise we find a quasispecies). Furthermore, Wahl et al. show that if susceptible host cells coexist with a resident lysogen-phage population with some lysogeny propensity ϕ¯r, a phage-lysogen pair of a variant with higher lysogeny propensity ϕ¯i>ϕ¯r can always invade on this equilibrium. Hence, in this equilibrium, the dominant phage will be the one with the highest equilibrium lysogeny propensity, ϕ¯=1 .

As has been demonstrated previously (Stewart and Levin, 1984; Wahl et al., 2018), the ‘coexistence equilibrium’ is stable only if phages and lysogens can invade on a susceptible population at carrying capacity (i.e. the condition in Equation A.16 is violated), and susceptible cells can invade on the phage-lysogen population in equilibrium (i.e. the condition in Equation A.20 is violated). Hence, susceptible cells, lysogens and phages all coexist precisely if

(A.21) (1-α)(B-1)a<δ<(B-1)a.

If the effective burst size B>1 (a necessary condition for the phage to be viable), this can be rewritten as

(A.22) (1-α)<δ(B-1)a<1.

Because the spontaneous induction rate of lysogens, α, is small (α<0.01, see Table 1), the condition in Equation A.22 is very specific. Susceptible cells, lysogens and phages coexist only if the exponential growth rate of a lytic phage spreading in a susceptible population at carrying capacity, (B-1)a-δ, is positive but very small, that is, if the epidemic is viable but only barely so. In reality, however, most phage epidemics are characterised by a high infectivity, mainly because of a large burst size (De Paepe and Taddei, 2006). Therefore, the condition in Equation A.22 is rarely satisfied, and for most phages we should instead expect to converge to equilibria of the third type (S¯=0,L¯>0,P¯>0).

This observation has consequences for the selection pressures on phage variants over the course of a typical epidemic. As soon as the pool of susceptible host cells is depleted, competition between the different phage variants vanishes and the relative frequency of the variants freezes (see Figure 2—figure supplement 1 for an illustration). Under these conditions, no infections take place and hence there is no selection on the lysis-lysogeny decision. We conclude that the evolution of the lysis-lysogeny decision of typical phages requires regular perturbations away from equilibrium conditions.

Appendix 3

Derivation of the evolutionarily stable strategy (ESS) under the serial-passaging regime

In the main text, we present simulation results of the evolution of the lysis-lysogeny decision of phages exposed to a serial-passaging regime, both when arbitrium communication is excluded from the model (Figure 2), and when it is included (Figures 3 and 4). These simulations show that, after many passaging episodes, a single variant dominates the phage population (accompanied by its quasispecies, see Figure 2C, Figure 3B). We therefore assume that, within the parameter range of interest, a pure Evolutionarily Stable Strategy (ESS) exists. In this section, we derive analytical expressions for the ESS. We first provide a definition of the ESS under a serial-passaging regime, and give a general description of how this ESS can be found (section A3.1). Then, we apply this general approach to derive the ESS of phages that differ in a constant lysogeny propensity, ϕ (absence of arbitrium communication, section A3.2), and the ESS of communicating phages that differ in their response threshold θ (section A3.3).

A3.1 General approach

Consider a population of phages under a serial-passaging regime with long time T between passages. At the start of each passaging episode, a fraction D of the phages is taken from the end of the previous episode, and is added to a ‘fresh’ population of susceptible bacteria at carrying capacity. This procedure is repeated over many episodes. Within each episode, the dynamics of the susceptible bacteria, lysogens, phages, and arbitrium are described by Equation A.7–A.12.

An evolutionarily stable strategy (ESS) is defined as a strategy that cannot be invaded by any other strategy that is initially rare. To find the ESS, we therefore consider a scenario where an invader phage variant attempts to invade a resident phage variant. Below, we specify what it means for a phage variant to be able to invade in a resident phage population under the imposed serial-passaging regime.

Envision a resident population consisting of an isogenic phage population that has gone through many passaging episodes. Over time, the dynamics within these episodes have converged to a repeatable trajectory characterised by Pr(t), the resident phage density over time, S(t), the density of susceptible bacteria over time, and Lr(t), the density of lysogens over time. At the start of one episode, now suddenly introduce a second phage with its own (possibly different) strategy. Crucially, the initial density of this invading phage Pi(0) is infinitesimally small. Consequently, during the first passaging episode the dynamics of the resident phage and the bacteria (Pr(t), S(t), and Lr(t)) are not affected by the invader. The invader is able to invade precisely if at the end of the first episode its frequency has increased relative to that of the resident, that is, if Pi(T)/Pr(T)>Pi(0)/Pr(0).

Note that the relative frequency of an invader with exactly the same strategy as the resident does not change during an episode. Suppose that such an invader is the best-performing invader under the environment set up by the resident; then this implies that no invader can increase in frequency over a passaging episode, and therefore the resident strategy must be an ESS. In other words, in a resident phage population consisting of the ESS only, the ESS itself is the optimal strategy for an invading phage variant, that is, the ESS is the optimal response to itself.

What does it take to be the best-performing invader? To answer this question, we consider the dynamics within a single passaging episode in more detail.

If the time between passages T is long (see below for exact condition), and the parameter conditions are such that the system converges to an equilibrium with S=0,P>0,L>0 (typical parameter values, see Appendix 2), the dynamics within an episode can be divided in three distinct phases (see Figure 2—figure supplement 1):

  1. Epidemic phase. A substantial population of susceptible bacteria (S>0) supports an ongoing epidemic. Free phages and lysogens are formed through infection of susceptible bacteria.

  2. Transition phase. The population of susceptible cells has collapsed (S0). The lysogen population expands to fill up the space left behind by lysed cells. Free phages particles can no longer infect susceptible cells, and disappear through decay and adsorption to lysogens.

  3. Equilibrium phase. The composition of the population is well-characterised by an equilibrium of ‘type 3’ (see Equation A.17). There is a small but consistent population of free phages that originates from lysogens through spontaneous induction. The distribution of phage variants in the free phage population is a direct representation of the relative frequency of the variants in the lysogen population (Equation A.18).

Let TE be the time at which the susceptible population collapses, that is, the end of the epidemic phase. If the time between passages T is sufficiently larger than TE, the passage takes place during the equilibrium phase. The relative frequency of phage variants in the passaged sample then directly reflects the relative frequency of the corresponding lysogens. Since lysogens are only differentially formed through infection dynamics (and not through lysogen replication, which happens at the same rate for all lysogen variants), the relative frequency of the different lysogens is established during the epidemic phase and does not change afterwards. The performance of an invading phage can hence be measured by its lysogen production between t=0 and t=TE.

The dynamics of S(t), and consequently TE, are determined by the resident phage: highly virulent resident phages (that cause many lytic infections, for instance because of a low ϕ-value) deplete the susceptible cell population faster than less virulent residents. The optimal invader under a certain resident is the phage variant that produces the most lysogens during the limited window of opportunity that it is offered by the environment set up by the resident. Since the ESS is the optimal response to itself, it is the strategy that, as an invader, produces the most lysogens during an epidemic phase set up by itself. We will use this reasoning to identify the ESS.

A3.2 Evolutionarily stable lysogeny propensity of non-communicating phages

Consider the restricted model in which arbitrium communication is excluded and phages are characterised by a fixed lysogeny propensity ϕ. To find the ESS under this scenario, we take the following steps:

  1. Derive how the duration of the epidemic, TE, depends on the ϕ-value of a resident phage population.

  2. Find the optimal ϕ given a fixed value of TE, that is, the ϕ-value that yields a maximal lysogen density at time TE.

  3. Combine 1. and 2. to find the ESS: the ϕ-value ϕ* that maximises its lysogen density at time TE(ϕ*), the duration of the epidemic as dictated by its own ϕ-value, ϕ*.

This approach is summarised in Box 1. Below, we provide the full analysis.

A3.2.1 Simplifying assumptions

To make the model analytically tractable, we make the following simplifying assumptions (based on the typical infection dynamics, see Figure 2A):

  1. Bacterial growth, decay of free phages, and induction of lysogens are considered to be much slower than the phage infection dynamics. We hence ignore these processes when describing the epidemic phase.

  2. The epidemic ends when all susceptible cells have been infected. In other words, we solve TE from t=0TEBaSPr dt=1 (where this one represents the carrying capacity in the non-dimensionalised units).

  3. The density of lysogens during the epidemic is small, hence NS.

  4. The susceptible population remains relatively constant for some time, after which it rapidly collapses. We approximate these dynamics with a block function, setting S=1 for tTE, and S=0 for t>TE.

Under these assumptions, the dynamics of the resident phage population for the period 0tTE are described by:

(A.23) dPrdt=(Ba(1-ϕr)S-aS)Pr=(B(1-ϕr)-1)aPr,

with solution

(A.24) Pr(t)=Pr,0e(B(1-ϕr)-1)at=Pr,0e(η-ϕr)Bat,

where Pr,0Pr(0) is the initial density of resident phages and we have introduced η:=1-B-1. Note that the description of the infection dynamics in Equation A.23 is meaningful only if early in the epidemic the phage population indeed grows exponentially, that is, if ϕr<η. For default parameter settings, this upper bound on ϕr is well above the ϕ-values that are typically selected (ϕ0.04 [Figure 2C], while η=0.5 [Table 1]), indicating that this assumption is reasonable.

A3.2.2 Duration of the epidemic given a resident phage

First, we derive how the duration of the epidemic, TE, depends on the lysogeny propensity of a resident phage variant ϕr. To find TE(ϕr), we substitute Equation A.24 into assumption 2:

t=0TEBaSPr(t) dt=t=0TEBaPr,0e(ηϕr)Bat dt=Pr,0(ηϕr)(e(ηϕr)BaTE1),

and then equate this integral to one to find

(A.25) TE(ϕr)=1Ba(η-ϕr)log(η-ϕrPr,0+1).

Note that the density of the resident phage at time TE is now given by

(A.26) Pr(TE)=Pr,0e(η-ϕr)BaTE=Pr,0+η-ϕr.

Therefore, the expression for TE (Equation A.25) can also be read as:

(A.27) TE(ϕr)=1Ba(η-ϕr)log(Pr(TE)Pr,0).

Equation A.27 will prove useful later for the derivation of the ESS.

A3.2.3 Optimal invader strategy given a fixed duration of the epidemic

Next, we ask what invader lysogen propensity ϕi,opt maximises the invader’s lysogen production, Li(TE), if the duration of the epidemic TE is fixed. For 0tTE, the dynamics of Li(t) are described by

(A.28) dLidt=ϕiBaPi.

Since Pi(t)=Pi,0e(η-ϕi)Bat (see Equation A.24) and Li(0)=0, we can now solve

(A.29) Li(t)=Pi,0ϕiη-ϕi(e(η-ϕi)Bat-1).

To find the ϕi-value that maximises Li(TE), we take the derivative of Equation A.29 with respect to ϕi:

(A.30) Li(TE)ϕi=Pi,0(η(e(η-ϕi)BaTE-1)(η-ϕi)2-ϕiBaTEe(η-ϕi)BaTEη-ϕi).

To find ϕi,opt, we should hence solve

(A.31) Pi,0η-ϕi(η(e(η-ϕi)BaTE-1)η-ϕi-ϕiBaTEe(η-ϕi)BaTE)=0.

Unfortunately, Equation A.31 cannot be solved analytically. We can however simplify Equation A.31 by noting that for sufficiently small ϕi, (η-ϕi) is of order 0.1-1, while TE is typically of order 1, and Ba is of order 10-1000 (Table 1). Hence, e(η-ϕi)BaTE is typically 1, and Equation A.31 can be approximated by

(A.32) Pi,0e(η-ϕi)BaTEη-ϕi(ηη-ϕi-BaϕiTE)=0.

From Equation A.32 we find

(A.33) ηη-ϕi=BaϕiTEBaTEϕ2-ηBaTEϕi+η=0ϕi,opt(TE)=12η±12η2-4ηBaTE.

A3.2.4 The ESS

Equation A.25 and its alternative formulation Equation A.27 describe how the duration of the epidemic TE depends on the lysogeny propensity ϕr of the current resident, while Equation A.33 gives the value ϕi,opt that maximises the lysogen production of an invader during an epidemic of a fixed duration TE. The ESS is now given by the value ϕ* that is ‘optimal’ as defined by Equation A.33, when it itself is the resident and hence dictates TE(ϕ*). Combining Equation A.33 and Equation A.27 we find

ϕ*=12η±12η2-4η(η-ϕ*)log(Pr(TE)/Pr,0)(ϕ*-12η)2=14η2-η(η-ϕ*)log(Pr(TE)/Pr,0)(ϕ*)2-ϕ*(η+ηlog(P(TE)/P0))+η2log(Pr(TE)/Pr,0)=0,

from which we can solve:

(A.34) ϕ*=12η(1+1log(Pr(TE)/Pr,0))±12η(1-1log(Pr(TE)/Pr,0)).

Of these two solutions, ϕ+=η is an asymptote at which our approximation no longer holds (remember that we previously demanded that ϕr<η to ensure initial spread of the infection). Hence, ϕ* should be given by ϕ-:

(A.35) ϕ*=ηlog(Pr(TE)/Pr,0).

Although Equation A.35 seems to provide an elegant equation for the ESS, it still depends on Pr(TE) and Pr,0. If the interval between passages is sufficiently long, the phage density at the end of a passaging episode will be given by Equation A.17 and hence

(A.36) Pr,0=Dα(1-α)δ+a(1-α),

where D is the dilution factor of phages upon passaging. While the value of Pr,0 does not depend on the lysogeny propensity ϕr, Pr(TE) does (see Equation A.26). Substituting Pr(TE)=Pr,0+η-ϕ* yields

(A.37) ϕ*=ηlog(1+η-ϕ*Pr,0).

This equation cannot be solved analytically. However, we can make a reasonable approximation of Equation A.37 by considering the differences in orders of magnitude of the terms within the logarithm. As argued above, (η-ϕ*) is generally of order 0.1-1, while typical values of Pr,0 are several orders of magnitude smaller (Pr,010-5). Therefore, we can approximate the logarithm in Equation A.37 by

log(1+η-ϕ*Pr,0)log(1Pr,0).

Using this approximation, we find

(A.38) ϕ*=ηlog(1Pr,0),

which is also presented in the main text (Equation 5). Equation A.38 and Equation A.36 were used to find the analytical predictions shown in Figure 2D.

A3.3 Evolutionarily stable response threshold of communicating phages

Next, consider a population of phages that do engage in arbitrium communication, again under a serial-passaging regime with sufficiently long time between the passages. Below, we use an approach similar to section A3.2, but more general, to derive the evolutionarily stable arbitrium response threshold, θ*. We take the following steps:

  1. Describe the dynamics of an invading phage and its corresponding lysogens in an environment dictated by a resident phage.

  2. Find the optimal invader response threshold under a fixed resident response threshold, that is, find the θ-value that maximises the invader’s lysogen production at time TE when the dynamics of susceptible cells (and hence TE) are determined by a fixed resident phage.

  3. Determine the ESS, θ*, as the optimal response to itself: the optimal invader response threshold (as found in step 2) if that same response threshold is the resident strategy.

We found that the results below are best understood in terms of the non-scaled model; in particular the (non-scaled) burst size of the phages turns out to be an important parameter. Therefore, the derivations below are presented for the dimensionalised equations Equation A.1–A.6.

A3.3.1 Simplifying assumptions

To make the model tractable, we again make a few simplifying assumptions:

  1. As in section A3.2, we assume that there is a separation of time scales between the infection dynamics of the phages and the reproduction of the bacteria, spontaneous phage decay and lysogen induction. Hence, when describing the epidemic phase we ignore these other processes.

  2. Additionally, we assume that there also is a separation of time scales between the production of arbitrium through infections (first term in Equation A.4) and its uptake and degradation by cells (second term in Equation A.4). We ignore the uptake and degradation of arbitrium during the early epidemic, such that the increasing arbitrium concentration reflects the decrease of the susceptible cell density because of infections.

  3. We assume that communicating phages switch from a completely lytic strategy (φ(A)=0) to a completely lysogenic strategy (φ(A)=1) once the arbitrium concentration exceeds the phages’ response threshold. This is in line with observations from simulations, where we find that phage variants with ϕmax=1 dominate the population for a wide range of parameter values (Figure 4A).

The assumptions above are less strict then the assumptions made in section A3.2. In particular, we no longer assume that the density of susceptible cells, S(t), remains constant for the duration of the epidemic 0tTE. Rather, for the derivations below it suffices to assume that S(t) is a declining function which is completely determined by the resident phage, and that S(t) is sufficiently close to 0 after the epidemic, that is, for times t>TE (where TE still depends on the characteristics of the resident phage).

It will be useful to refer to the time when the resident and invader switch to lysogeny as τr and τi, respectively. Given particular dynamics of the susceptible cell density and the arbitrium concentration dictated by a resident phage population, there is a direct relation between τi and θi (the arbitrium response threshold of the invader). Keep in mind, however, that this relation changes if the resident phage is changed.

A3.3.2 Dynamics of the invading phage and its corresponding lysogens

First, we describe how the dynamics of an invading phage variant depend on the resident phage population. Remember that we consider an invader phage variant that starts off at infinitesimally small density, and attempts to invade an isogenic resident phage population that has already converged to a repeatable trajectory of Pr(t), Lr(t), S(t), and A(t) per passaging episode. Under these conditions, the dynamics of S(t), N(t)=S(t)+Lr(t) and A(t) over the first passaging episode do not depend on the switch time τi of the invader, but only on the switch time of the resident phage, τr. Based on the assumptions formulated above, the ODEs for the density of the invading phage and its corresponding lysogens can be written as

(A.39)dPi(t|τi,τr)dt={BbaS(t|τr)Pi(t|τi,τr)aN(t|τr)Pi(t|τi,τr),(t<τi)aN(t|τr)Pi(t|τi,τr),(tτi)(A.40)dLi(t|τi,τr)dt={0,(t<τi)baS(t|τr)Pi(t|τi,τr),(tτi)

where vertical lines are used to indicate which phage charactistics the trajectories of variables depend upon. Equation A.39–A.40 capture the switch from a completely lytic infection strategy (for tτi), in which new phage particles are produced through infection but no lysogens are formed, to a completely lysogenic strategy (for t>τi), in which no new phage particles are produced but all infections result in the production of lysogens (see assumption 3). Remember that we here use the dimensionalised equations, so B is the burst size, a the rate of adsorption of phages to bacterial cells (irrespective of whether they are susceptible or lysogen), and b the probability that adsorption to a susceptible cell leads to an infection.

The solution to Equation A.39 can be written as

(A.41) Pi(t|τi,τr)=Pi,0×{exp(Bba0tS(t|τr)𝑑t-a0tN(t|τr)𝑑t),(t<τi)exp(Bba0τiS(t|τr)𝑑t-a0tN(t|τr)𝑑t),(tτi)

as is easily verified by differentiating this solution with respect to t.

As before, the performance of the invading phage variant is determined by its lysogen production during the epidemic phase, that is, between t=0 and t=TE. At any time t, the density of invader lysogens is

(A.42) Li(t|τi,τr)={0,(t<τi)τitbaS(t|τr)Pi(t|τi,τr)dt.(tτi)

Invading phage variants are selected on their lysogen density at the end of the epidemic, Li(TE|τi,τr). Once the epidemic phase has ended (tTE), no new phage particles or lysogens are formed through infection. Hence, any reasonable switching time must obey τi<TE. Furthermore, since S(t)0 for any time tTE,

(A.43) Li(TE|τi,τr)=τiTEbaS(t|τr)Pi(t|τi,τr)𝑑t=τibaS(t|τr)Pi(t|τi,τr)𝑑t

which we will denote Λi(τi,τr).

A3.3.3 Optimal invader strategy given some resident phage

The optimal invader strategy τi given S(t|τr) and N(t|τr) is the one that maximises Λi(τi,τr). To find this optimal strategy, we differentiate Equation A.43 with respect to τi:

(A.44) Λi(τi,τr)τi=-baS(τi|τr)Pi(τi|τi,τr)+τibaS(t|τr)Pi(t|τi,τr)τi𝑑t.

The derivative in the integrand can be calculated from Equation A.41 (noting that, inside the integral, tτi):

(A.45) Pi(t|τi,τr)τi=τi(Bba0τiS(t|τr)dt)Pi,0exp(Bba0τiS(t|τr)dta0tN(t|τr)dt)=BbaS(τi|τr)Pi(t|τi,τr).

Inserting the last expression into Equation A.44 yields

(A.46) Λi(τi,τr)τi=baS(τi|τr)Pi(τi|τi,τr)+Bb2a2S(τi|τr)τiS(t|τr)Pi(t|τi,τr)dt=baS(τi|τr)Pi(τi|τi,τr)+BbaS(τi|τr)Λi(τi,τr).

The terms in Equation A.46 have a clear interpretation. By taking the derivative of Λi(τi,τr) to τi, we are implicitly comparing one possible invading phage variant (phage 1) that switches at time t=τi to a second invading phage variant (phage 2) that switches ever so slightly later, at t=τi+dτ. Equation A.46 says that the lysogen density of these two variants at the end of the epidemic will differ because of two effects: On the one hand (first term) phage 2 will have a lower lysogen density than phage 1 because it does not produce lysogens in the time interval from τi to τi+dτ. The damage is -baS(τi|τr)Pi(τi|τi,τr)dτ lysogens per volume. On the other hand, phage 2 will have a higher higher lysogen density because it produces additional free phages in the time interval from τ to τ+dτ, which results in additional lysogens in the rest of the epidemic. As a result, throughout the rest of the epidemic the second phage has (1+BbaS(τi|τr)dτ) times as many phages as the first phage variant, and therefore produces an additional number of BbaS(τi|τr)Li(τi,τr)dτ lysogens per volume.

The optimal invading phage variant given a resident phage is the variant with the value of τi,opt(τr) for which the two terms in Equation A.46 cancel precisely:

(A.47) Pi(τi,opt(τr)|τi,opt(τr),τr)=BΛi(τi,opt(τr),τr).

That is, the optimal invader switches precisely when its phage density is equal to its total eventual lysogen production multiplied by the burst size B.

We may rewrite Equation A.47 as

(A.48) BEi(τi,opt(τr)|τr)=1,

where Ei(τi|τr)Λi(τi,τr)/Pi(τi|τi,τr) is the number of lysogens eventually produced per phage of the invader phage variant present at time τi. Ei(τi|τr) can be interpreted as a kind of ‘exchange rate’, expressing the value of a single phage at time τi in the currency of lysogens. This suggests another way of phrasing the results above, where we compared two phage variants of which phage 2 switched slightly later than phage 1: During the time interval from τi,1 to τi,2=τi,1+dτ, both competing invading phage variants infect baS(τi,1)Pi(τi,1|τi,1,τr)dτ susceptible bacteria per volume. Phage 1 directly converts these infected bacteria into lysogens. Phage 2 instead converts each of them into B additional phages. Whether this is a good idea depends precisely on whether increasing the phage density by B phages per volume will, during the rest of the epidemic, result in an increased lysogen density of more than 1 lysogen per volume. That is, phage 2 is the better invader precisely if BEi(τi|τr)>1, while phage 1 is the better invader if BEi(τi|τr)<1. Again we see that the optimal invader must obey Equation A.47 and Equation A.48.

A3.3.4 The ESS

To find the ESS, we ask what phage variant is the optimal response to itself, that is, what phage variant satisfies

(A.49) τi,opt(τ*)=τ*.

In other words, the ESS must obey a special case of Equation A.47 and Equation A.48:

(A.50) Pi(τ*|τ*,τ*)=BΛi(τ*,τ*) or BEi(τ*|τ*)=1.

Importantly, in this case the resident and invader behave identically, so that the exchange rate of the resident must be the same as that of the invader. That is,

(A.51) Ei(τ*|τ*)=Er(τ*|τ*)=Λr(τ*)/Pr(τ*|τ*),

where Λr(τr) is the total density of lysogens eventually produced by a resident with switching time τr. Combined with Equation A.50 this results in:

(A.52) BΛr(τ*)/Pr(τ*|τ*)=1.

Hence, the ESS is the strategy that, when it is the only phage variant present, switches from the lytic to the lysogenic cycle precisely when the density of free phage particles it has is equal to the burst size times the density of lysogens it will still produce in the remainder of the active epidemic (see explanation in the previous section).

So far, we have expressed the results for the ESS as a switching time τ*. In reality, however, the communicating phages switch when a certain threshold arbitrium concentration θ* is reached. As a last step, we therefore have to relate the terms in Equation A.52 to the arbitrium concentration. Under our simplifying assumptions, the arbitrium dynamics between time t=0 and t=τr are described by

(A.53) A(t|τr)=cbaS(t|τr)Pr(t|τr),

where c is increase in arbitrium concentration per infection. The total arbitrium concentration at time t is hence given by

(A.54) A(t|τr)=0tcbaS(t|τr)Pr(t|τr)𝑑t,

which can be written as A(t|τr)=cIr(t|τr), where

(A.55) Ir(t|τr)=0tbaS(t|τr)Pr(t|τr)𝑑t

is the infection density: the number of infections that has occurred per volume at time t.

To express the ESS in terms of the arbitrium concentration, we first show that Pr(τ*|τ*) is approximately proportional to Ir(τ*|τ*). In general, the resident phage density obeys an equation equivalent to Equation A.39 (even though this equation was originally written down for the invading phage). For the time period t<τr, the solution of this equation can be expressed as

(A.56) Pr(t|τr)=Pr,0+0t(Bb1)aS(t|τr)Pr(t|τr)dt=Pr,0+(Bb1)Ir(t|τr)

Provided that the initial phage density Pr,0 is negligible compared to the phage density at time τ*, we find that

(A.57) Pr(τ*|τ*)(B-b-1)Ir(τ*|τ*).

Next, we use that the epidemic will eventually consume (almost) all susceptible bacteria. (Note that this is equivalent with our earlier assumption that S(t)0 for t>TE.) Hence, we must have that

(A.58) Λr(τ*)S(0)-Ir(τ*|τ*).

If we insert Equation A.57 and Equation A.58 into Equation A.47 and solve for Ir(τ*|τ*), we arrive at

(A.59) Ir(τ*|τ*)=S(0)2-(bB)-1.

That is, the ESS switches when the infection density obeys Equation A.59. This implies that the ESS should have the threshold

(A.60) θ*=cIr(τ*|τ*)=cK2-(bB)-1,

where we have substituted S(0)=K, the carrying capacity of the bacteria. Equation A.60 is also presented in the main text (Equation 6). This equation was used to provide the analytical estimates shown in Figure 4B.

Data availability

All data were obtained through computer simulation. Scripts to run these simulations, simulated data, and analysis scripts are available at GitHub: https://github.com/hiljedoekes/PhageCom (copy archived at https://archive.softwareheritage.org/swh:1:rev:8124bcd18e18dd03b94a315b7694af3ed2e4a002/).

References

    1. Hershey AD
    (1946)
    Mutation of bacteriophage with respect to type of plaque
    Genetics 31:620–640.

Decision letter

  1. Samuel L Díaz-Muñoz
    Reviewing Editor; University of California, Davis, United States
  2. Aleksandra M Walczak
    Senior Editor; École Normale Supérieure, France
  3. Samuel L Díaz-Muñoz
    Reviewer; University of California, Davis, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

The recent discovery of phage communication allowing collective decision making (Erez et al., 2017), raised many questions about the social lives of viruses. This paper uses a mathematical modeling approach to determine the conditions under which a small molecule communication system among phages governing the lysis-lysogen decision would evolve. The paper, the first to employ theory on the topic, demonstrates the evolutionary advantages of phage communication over a bet-hedging strategy, which was previously believed to govern life cycle decisions in viruses.

Decision letter after peer review:

Thank you for submitting your article "Repeated outbreaks drive the evolution of bacteriophage communication" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Samuel L Díaz-Muñoz as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Aleksandra Walczak as the Senior Editor.

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

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, when editors judge that a submitted work as a whole belongs in eLife but that some conclusions require a modest amount of additional new data, as they do with your paper, we are asking that the manuscript be revised to either limit claims to those supported by data in hand, or to explicitly state that the relevant conclusions require additional supporting data.

Our expectation is that the authors will eventually carry out the additional experiments and report on how they affect the relevant conclusions either in a preprint on bioRxiv or medRxiv, or if appropriate, as a Research Advance in eLife, either of which would be linked to the original paper.

Summary:

This paper uses a modeling approach to determine the conditions under which a small peptide-mediated communication as a strategy informing lysis-lysogeny decision would arise in viruses. This topic arises from the recent discovery of phage communication by Erez et al., 2017. It starts using traditional differential equation system to investigate the lysis-lysogeny decision in a bet hedging context. However, most of the manuscript goes on to describe a modeling strategy mimicking serial passaging to examine the bet hedging strategy and the communication strategies, including examining how they work in a head to head competition (when does the communication strategy invade). Finally, the paper examines the threshold of proportion of cells infected at which phages are predicted to switch to lysogeny, coming up with 0.5 as a consistent number across a wide parameter swath. The manuscript demonstrates how such a communication mechanism can be fitter than a bet-hedging strategy where some fixed fraction of infections result in lysogeny.

In sum, this manuscript is very clear, well-argued, well-written, and it was a pleasure to read it. This paper is a brilliant contribution not only for being (surprisingly) the first to address this topic from a theoretical perspective, but for its thoughtfulness, rigor, completeness, and the connection of the theoretical and the biological. There are a few areas for improvement, but this manuscript should be an outstanding contribution, certainly to the emergent sociovirology literature, but also to a wide audience including microbiologists, virologists, ecologists, evolutionary biologists, and behaviorists.

Please pay very special attention to the “Revisions expected in follow-up work” section below in revising your manuscript, the other revisions are there for our consideration and are given with the intention of improving the manuscript.

Revisions for this paper:

1) Some assumptions of the model should be better justified and/or discussed, for the reader to assess the generality of the conclusions obtained in the manuscript. In particular:

a) Subsection “Model” paragraph three: It would be good to comment on how realistic the stepwise assumption is for phi(A).

b) Results third paragraph: This discussion strongly relies on the fact that cells stop dividing as carrying capacity is reached. In reality, cells may die at a certain rate, or some cells may leave the system or be lost, implying that some rate of division could be sustained at equilibrium. It would be really good to discuss this possibility and the impact it may have on results. The serial passaging scenario considered afterwards is a special way of including such an effect, but at discrete and periodic times.

Fourth paragraph: It would be important to further motivate the serial passaging scenario that is chosen. Is it related to actual or potentially realistic experiments? Are the conclusions obtained robust to modifications of this scenario?

2) Results second paragraph: I completely agree with the statements and rationale outlined in this paragraph. However, the scenario described in paragraph four should then be part of the Materials and methods. I still think it can be included in the Results as is, but needs to be included in the Materials and methods.

3) Some elements of Discussion should be specified:

(3a) Paragraph four: Two mutations occurring in a very short lapse of time seems rather unlikely a priori. Could the authors comment on whether scenarios where one happens and then the other would be favorable? Related to mutations: In the passaging is there selection with regard to the phage mutations or are they likewise proportionally represented in the next passage? If they are rare in the system, they could be lost in the dilution. How would this affect the potential evolution and impact the model outcomes?

b) Could there be other scenarios under which the production of arbitrium by lysogens is more useful?

c) Final paragraph: It would be great if the authors could give some indications about how the current model could be adapted to these other cases and what insight it may bring. This would add value to the manuscript by potentially making its scope broader.

4) The readability of Figure 1A, which I think needs revision. I can see the authors likely spent considerable effort on it, but it remains very difficult to parse. Is it really necessary to have four different kinds of arrows? As the eye wanders around the diagram, it's not clear where to start or where to end. Perhaps it would help to organizing cycles more neatly into circles? The many virus pictograms add visual clutter without aiding clarity. I don't have a clear recipe, but I'd encourage the authors to solicit input from colleagues outside their field so as to make this figure maximally accessible, especially for a broad-audience-journal like eLife. Elsewhere the authors do a truly admirable job, reminding the reader what their parameters are etc., so it would be a shame for the reader to stumble at Figure 1A.

5) The discussion is brief, but excellent. It really covers a lot of the questions that I wanted to know before reading the paper and does an amazing job of connecting model results to biology and making predictions to test empirically.

Revisions expected in follow-up work:

1) My major "science concern" was spurred right from the start by the sentence in the Abstract, "our model predicts the selection of phages that switch infection strategy when half of the available susceptible cells have been infected". Seeing as it is notoriously difficult for a virus to estimate what fraction of cells were infected, this raised a worry that persisted through reading the entire paper – until the very last paragraph, where it was finally discussed. My suggestions would be to do this (a) sooner, and (b) better.

a) Sooner: I'd recommend explicitly commenting on the issue already in the Introduction and at least promise that it will be discussed.

b) Better: “The arbitrium concentration during the early epidemic then is a direct reflection of the fraction of susceptible cells that have so far been infected”

The readout is a direct reflection only of how many cells were lysed. In the model, yes, this is the same as fraction. But as the authors acknowledge at the very end this requires the assumption that bacterial cells were at carrying capacity K. In natural conditions, the carrying capacity can differ due to environmental factors – "unknown" to the virus.

Far from being a minor worry, the issue is directly relevant for the authors' comparison between the communicating and the bet-hedging strategy. Under the communicating strategy, if the carrying capacity were modulated to e.g. half of its usual value, it seems that the virus would never switch to the lysogenic phase – a crippling blow to fitness. In contrast, the bet-hedging strategy seems vastly more robust to such a modulation. This suggests that the conclusions could be changed quite strongly if the carrying capacity were picked randomly (within some range) between passaging trials.

It sounds like there is a prediction here.

I do not intend this comment as negative / undermining authors' findings. Quite the opposite, I think by adding this analysis the authors could strengthen their argument, demonstrating the utility of a quantitative framework like theirs, and perhaps even make additional predictions. Beyond what level of (passage-to-passage) variability of K the bet-hedging strategy becomes favored? Addressing this comprehensively would be a project in itself; but generating a figure panel illustrating this in simulations would strengthen the authors' case for the benefits of modeling work.

An optional point:

Since the cell-mediated uptake and degradation of arbitrium is also mediated by bacterial cells, and is thus also dependent on bacterial cell density, it's possible that this might mitigate the problem (making the arbitrium concentration less strongly dependent on K than naïve linearity). I don't know if it's true, but if it is, it would be cool, and worth highlighting. Unless I'm mistaken, this is not currently discussed. I think that could strengthen the story.

2) Figure 3: In the competition between bet-hedging and communicating variants, what is the impact of initial conditions? Moreover, here, only the 2 top variants of each type are competed against each other. Is this sufficient to conclude that communicating variants will always be selected by competition? If binary comparisons are sufficient, this should be stated and explained. If not, then it would be important to show explicitly what happens when various phages of both types compete together.

3) Scripts really should be openly available unless there is a compelling reason not to do so, especially in this case where the results depend so closely on the code. I would strongly urge the authors to instead publish all simulation scripts alongside the manuscript, both for archiving purposes but also to facilitate access.

4) The parameters tested should be in the main text. Note that, other than the construction of the equations according to life history parameters, this is the only connection to "biology" that non-theoretical researchers will have to the manuscript. Table A1 should be in the main text. The table should list the conventional units of the parameters. For instance does burst size of 2 refer to 2 virions? This value would normally be in the 10's-100's for most phages. Some of these parameters don't seem reasonable at first glance. For instance, the adsorption rate is usually something on the order of 10-10 to 10-12. I suspect this could be because these are "scaled" values. If the scaled values are made for model convenience they should also be listed separately. This way the connection between the model and the biology will be more evident to readers.

5) The idea of having to go away from a steady state for this to work is very appealing. It is also likely much closer to biological reality in this case. Other than models that specifically address serial passaging, is the approach of disturbing ESS's regularly a new concept? Is there a generalizable framework for this?

6) Some conclusions should be specified or clarified:

a) Figure 2 B-C: Some diversity seems to survive at steady state for T>5h. I believe this is what the authors refer to as "quasi-species" in the manuscript. It would be important to discuss and interpret this surviving diversity, and to explain whether it is going to survive forever or to slowly decay, and why.

b) Figure 4: Please explain the difference between the analytical prediction and the simulation results at small values of bB in Figure 4B. In particular, does any of the simplifying assumptions result in this discrepancy?

https://doi.org/10.7554/eLife.58410.sa1

Author response

Revisions for this paper:

1) Some assumptions of the model should be better justified and/or discussed, for the reader to assess the generality of the conclusions obtained in the manuscript. In particular:

a) Subsection “Model” paragraph three: It would be good to comment on how realistic the stepwise assumption is for phi(A).

To our knowledge, few details are currently known about the response function φ(A). Therefore, we had to make some simplifying assumption about the shape of this curve. We do agree that it is fair to acknowledge this, and to discuss this assumption. We now do so in the Discussion.

b) Results third paragraph: This discussion strongly relies on the fact that cells stop dividing as carrying capacity is reached. In reality, cells may die at a certain rate, or some cells may leave the system or be lost, implying that some rate of division could be sustained at equilibrium. It would be really good to discuss this possibility and the impact it may have on results. The serial passaging scenario considered afterwards is a special way of including such an effect, but at discrete and periodic times.

We think a misunderstanding has arisen here. The susceptible host cell population is depleted at equilibrium not because cells no longer divide, but because the susceptible host cells are outcompeted by lysogens that are immune to superinfection and hence do not experience the additional death rate induced by phage infections. We have altered the text to stress that the susceptible cells are outcompeted by lysogens.

The logistic terms used to describe bacterial growth in our model do not assume that cells stop dividing when carrying capacity is reached. Rather, they assume that at carrying capacity the division rate of cells is equal to the death (or loss) rate, such that a stable density of bacteria is maintained. In other words: the logistic growth rate is a net growth rate (division – death).

To clarify this point, let us show mathematically that an additional density-independent death rate can always be absorbed into a logistic term of net growth. Assume that the bacterial division rate is density dependent (and is hence described by a logistic-like term), while the bacterial death rate does not depend on density. The net growth of bacterial density B is then described by dBdt=dB(1bh)dB where b is the bacterial division (or birth) rate, h is the bacterial density at which the division rate is equal to zero, and d is the density-independent death rate. Now, note that we can rewrite this equation as: dBdt=(bd)BbB2h=(bd)B (1b(bd)h B) If we now define r = b/d and K = (1 – d/b)h, this equation is exactly equal to the logistic terms used in our model.

Fourth paragraph: It would be important to further motivate the serial passaging scenario that is chosen. Is it related to actual or potentially realistic experiments? Are the conclusions obtained robust to modifications of this scenario?

On the one hand, the serial passaging scenario was inspired by the set-up of actual serial passaging experiments (e.g., Bull et al., Evolution, 1993; Bollback et al., Mol Biol and Evol, 2007). On the other hand, it is also meant to mimic the large changes in the density of available susceptible host cells probably experienced by phages in real life situations. This is now stated more clearly in the Materials and methods.

In terms of robustness, several modifications to the serial passaging set-up can be imagined. We chose to study two: varying the carrying capacity of bacteria between cycles, and altering which entities are passaged. The first of these modifications was also motivated by another comment (“Revisions expected in follow-up work – 1b”) and is discussed in detail there. In summary, we found the results to be quite robust against these variations in bacterial carrying capacity (see new Figure 5 and subsection “Arbitrium communication is robust against variation in bacterial carrying capacity”). For the second modification, note that in our initial set-up only phages are transferred from one epidemic cycle to the next. To test the robustness of the results against alterations of this assumption, we repeated our simulations but now passaged a small sample of the full system (i.e., susceptible bacteria, lysogens, phages, and arbitrium). We found that the results were highly robust to this change in the set-up. For the case without arbitrium signalling these results are shown in the new (Figure 2—figure supplement 2 ), and for the case with arbitrium signalling they are shown in the new Figure 3—figure supplement 2.

2) Results second paragraph: I completely agree with the statements and rationale outlined in this paragraph. However, the scenario described in paragraph four should then be part of the Materials and methods. I still think it can be included in the Results as is, but needs to be included in the Materials and methods.

We understand this concern. We have moved the technical description of the serial passaging regime to the Materials and methods, while maintaining our explanation of why such a regime is necessary in the Results.

3) Some elements of Discussion should be specified:

(3a) Paragraph four: Two mutations occurring in a very short lapse of time seems rather unlikely a priori. Could the authors comment on whether scenarios where one happens and then the other would be favorable?

Indeed, it is not really necessary for these two mutations to take place within a very short lapse of time. A mutation that allows phages to escape superinfection immunity will always free up a new pool of susceptible cells, namely the lysogens that are now no longer protected against this altered phage. If this mutant phage variant however still responds to the arbitrium that was produced during previous infections of the original phage variant and this arbitrium is still present in a sufficiently high concentration, the phage likely makes the “wrong” lysis-lysogeny decision: to cause lysogenic infections and hence replicate slowly while it could replicate much more quickly through the lytic cycle. Under these circumstances, there will be selection pressure on the mutant phages to acquire additional mutations that alter the signalling peptide.

We have altered the discussion to describe this scenario.

Related to mutations: In the passaging is there selection with regard to the phage mutations or are they likewise proportionally represented in the next passage?

There is no selection during the passaging and phage variants are proportionally included in the passaged sample. This is described in the introduction of the serial passaging set-up (Materials and methods).

If they are rare in the system, they could be lost in the dilution. How would this affect the potential evolution and impact the model outcomes?

Since our model is a deterministic model of continuous densities and not a stochastic model of discrete individuals, loss by dilution is impossible. The density of particular strains may become very low but is never actually zero. This is a consequence of the choice of modelling framework, and obviously not a true representation of nature. If a certain phage variant is very rare in a natural system, it can of course be lost by dilution. However, since we consider phage variants that can be formed through mutations of other phages, such a loss could always be compensated by a mutation that reinstates the lost variant. How quickly such a mutant arises depends on population characteristics such as the effective population size. Demographic stochasticity (including loss by dilution) can hence quantitatively alter the evolutionary dynamics. However, as long as the dynamics are not completely dominated by such stochasticity (i.e., selection is not completely overshadowed by drift), the qualitative results such as the ESS will be the same.

b) Could there be other scenarios under which the production of arbitrium by lysogens is more useful?

This is a fair question, and we have given it considerable thought. The main reason that the production of arbitrium by lysogens is not beneficial in our modelling set-up is that lysogens and susceptible cells do not coexist for any significant amount of time. If there are no susceptible cells in the presence of lysogens, then there are no infection events and hence arbitrium produced by the lysogens cannot influence any lysis-lysogeny decisions. As we show in the analysis of the model’s equilibria and the parameter conditions for these equilibria to be stable, coexistence between susceptible cells and lysogens in equilibrium is very unlikely. However, coexistence could be attained under different circumstances. The most reasonable scenario is then that somehow there is a constant influx of susceptible cells into a population of lysogens.

We see two ways in which this could happen. Firstly, there might be an actual, physical inflow of susceptible cells. In natural systems, large-scale outbreaks could then correspond to phage colonization of a new area (e.g., a particular part of a human gut), while the low-rate constant inflow is caused by movement of cells in this area (e.g., bacterial migration within the gut).

Secondly, a low constant influx of susceptible cells into a population of lysogens could arise from the occurrence of loss-of-super immunity mutations in the prophages carried by the lysogens which render the lysogen vulnerable to the phage again. Since such mutations likely happen at low rates the “influx rate” of susceptible cells will be low, but such a process might allow a small population of susceptible cells to coexist with a population of lysogens. Since there will be only few of these susceptible cells at any given time, a lysogenic strategy is favoured when these cells are infected, and phages might hence benefit from arbitrium production by lysogens.

We now discuss both possibilities in the Discussion.

c) Final paragraph: It would be great if the authors could give some indications about how the current model could be adapted to these other cases and what insight it may bring. This would add value to the manuscript by potentially making its scope broader.

Thank you for this suggestion. We have added a short description of possible adaptations to our model and what we might learn from them.

4) The readability of Figure 1A, which I think needs revision. I can see the authors likely spent considerable effort on it, but it remains very difficult to parse. Is it really necessary to have four different kinds of arrows? As the eye wanders around the diagram, it's not clear where to start or where to end. Perhaps it would help to organizing cycles more neatly into circles? The many virus pictograms add visual clutter without aiding clarity. I don't have a clear recipe, but I'd encourage the authors to solicit input from colleagues outside their field so as to make this figure maximally accessible, especially for a broad-audience-journal like eLife. Elsewhere the authors do a truly admirable job, reminding the reader what their parameters are etc., so it would be a shame for the reader to stumble at Figure 1A.

We understand this comment and have tried to make the figure easier to digest. We removed many virus pictograms and restructured the diagram such that the figure can be read starting in the upper left corner. The infection event and subsequent lysis-lysogeny decision now take center stage. The lytic and lysogenic cycle are shown as two equally shaped ellipses, to stress the fact that these are two alternative life cycles. We have furthermore included a colour correspondence between the process influenced by the arbitrium concentration (the lysis-lysogeny decision) and the arbitrium itself.

Revisions expected in follow-up work:

1) My major "science concern" was spurred right from the start by the sentence in the Abstract, "our model predicts the selection of phages that switch infection strategy when half of the available susceptible cells have been infected". Seeing as it is notoriously difficult for a virus to estimate what fraction of cells were infected, this raised a worry that persisted through reading the entire paper – until the very last paragraph, where it was finally discussed. My suggestions would be to do this (a) sooner, and (b) better.

We understand this concern and thank the reviewer for his/her detailed consideration and suggestions (discussed further below).

a) Sooner: I'd recommend explicitly commenting on the issue already in the Introduction and at least promise that it will be discussed.

We agree that it is better to address this point earlier in the manuscript. We have now further investigated how accurate the information carried by the arbitrium signal needs to be for communication to evolve (see point 1b), and have added a sentence to the Introduction referring to these new results: “Finally, we investigate how reliable the arbitrium signal needs to be for such communication to evolve, and find that the results are remarkably robust against variation in the density of bacteria”.

b) Better: “The arbitrium concentration during the early epidemic then is a direct reflection of the fraction of susceptible cells that have so far been infected”

The readout is a direct reflection only of how many cells were lysed. In the model, yes, this is the same as fraction. But as the authors acknowledge at the very end this requires the assumption that bacterial cells were at carrying capacity K. In natural conditions, the carrying capacity can differ due to environmental factors – "unknown" to the virus.

Far from being a minor worry, the issue is directly relevant for the authors' comparison between the communicating and the bet-hedging strategy. Under the communicating strategy, if the carrying capacity were modulated to e.g. half of its usual value, it seems that the virus would never switch to the lysogenic phase – a crippling blow to fitness. In contrast, the bet-hedging strategy seems vastly more robust to such a modulation. This suggests that the conclusions could be changed quite strongly if the carrying capacity were picked randomly (within some range) between passaging trials.

It sounds like there is a prediction here.

I do not intend this comment as negative / undermining authors' findings. Quite the opposite, I think by adding this analysis the authors could strengthen their argument, demonstrating the utility of a quantitative framework like theirs, and perhaps even make additional predictions. Beyond what level of (passage-to-passage) variability of K the bet-hedging strategy becomes favored? Addressing this comprehensively would be a project in itself; but generating a figure panel illustrating this in simulations would strengthen the authors' case for the benefits of modeling work.

This is a great suggestion, and we have now done such follow-up work. In these new simulations, at the start of each passaging episode a random value of the carrying capacity is sampled from a γ distribution. Sampling from a γ distribution is convenient for 2 reasons: (1) it allows us to change the variance of the distribution while maintaining the same mean, and (2) values sampled from this distribution are always >= 0. We performed these simulations for a long between-passage time (T = 24 h) and for varying levels of carrying capacity noise (i.e., different values of the variance). The results of these new simulations are shown in a new Figure 5 and Figure 5—figure supplement 1.

Two things stand out when considering these new results. Firstly, a communication strategy with φmax = 1 is selected for simulations in which the carrying capacity varies with a coefficient of variation (CV = standard deviation / mean) up to 0.35. The result that phages are selected that switch from a fully lytic to a fully lysogenic strategy is hence surprisingly resistant to noise in the carrying capacity. We do see that the response threshold value that is selected slightly declines with increasing levels of noise. This makes sense: as also stated by the reviewer, any passaging episodes in which the carrying capacity is lower than the phages’ response threshold are disastrous for the virus because under these conditions no lysogens are produced. As the variation in the bacterial carrying capacity increases, the phages hence become “more prudent”: they switch strategies earlier in each outbreak to avoid not switching at all in some outbreaks.

Secondly, the results partly confirm the hypothesis put forward by the reviewer, but with an important nuance. At high levels of variation in the carrying capacity (CV > 0.35), phages are selected with φmax < 1. These phages hence show a bet-hedging strategy when the arbitrium concentration is above their response threshold value. Surprisingly, however, the response threshold value θ of the selected phages is never zero. We find this result for all levels of variation we tested, up to simulations in which the standard deviation on the carrying capacity is equal to its mean value (CV = 1). This suggests that even for very large variation in the bacterial carrying capacity a form of arbitrium communication is favoured over a completely bet-hedging strategy. Over the course of an outbreak, the selected phages still start with a lytic strategy, but do switch to a bet-hedging strategy at some low (but non-zero) arbitrium concentration.

To conclude, we find that the arbitrium communication system is more robust to variations in bacterial carrying capacity than we expected.

We included a short description of the set-up of these new simulations in the Materials and methods and discuss the results in a new section in the Results. Based on these new results, we also changed the corresponding paragraph in the Discussion.

An optional point:

Since the cell-mediated uptake and degradation of arbitrium is also mediated by bacterial cells, and is thus also dependent on bacterial cell density, it's possible that this might mitigate the problem (making the arbitrium concentration less strongly dependent on K than naïve linearity). I don't know if it's true, but if it is, it would be cool, and worth highlighting. Unless I'm mistaken, this is not currently discussed. I think that could strengthen the story.

This is an interesting idea. However, for the parameter conditions used here it unfortunately does not seem to hold true (illustrated in Author response image 1). This has two reasons. First, the epidemiological dynamics tend to be very fast, and hence the production of arbitrium (due to infection events) overrules its uptake and degradation early in the epidemic. Second, during the active epidemic the number of cells temporarily declines because susceptible cells get lysed and their place is only later taken over by lysogens. Hence, the total uptake rate of arbitrium during the crucial active epidemic phase is relatively low due to the low number of cells.

Although we like the idea, we decided not to include it in the manuscript because the effect does not seem to occur in the current model set-up.

Author response image 1
Short-term dynamics for lower values of the bacterial carrying capacity K.

2) Figure 3: In the competition between bet-hedging and communicating variants, what is the impact of initial conditions? Moreover, here, only the 2 top variants of each type are competed against each other. Is this sufficient to conclude that communicating variants will always be selected by competition? If binary comparisons are sufficient, this should be stated and explained. If not, then it would be important to show explicitly what happens when various phages of both types compete together.

The binary comparison shown in Figure 3C on its own is indeed not sufficient to conclude that communicating variants are always selected over bet-hedging variants. However, bet-hedging phage variants are also included in the large collection of possible variants in the simulations of Figure 3B. Strains with response threshold θι = 0 always have a lysogeny propensity of φmax,I, and hence do not respond to the arbitrium signal but rather employ a bet-hedging strategy. If any bet-hedging variant would be favoured over the communicating phages, it should hence dominate in Figure 3B. This does not happen, so from Figure 3B we can already conclude that the communicating strategy is favoured over the bet-hedging strategy. The one-on-one competition in Figure 3C was only added to underscore this conclusion.

To make clear that bet-hedging phages are included in the eco-evolutionary simulations with arbitrium communication, we now explain this in the Materials and methods. We have also rephrased part of the discussion of Figure 3B and 3C, to clarify that 3C is an example underscoring the conclusion already drawn from Figure 3B.

3) Scripts really should be openly available unless there is a compelling reason not to do so, especially in this case where the results depend so closely on the code. I would strongly urge the authors to instead publish all simulation scripts alongside the manuscript, both for archiving purposes but also to facilitate access.

We agree that publishing code is good practice, and apologise for initially taking the easy way out. We have now prepared the scripts for publication and made them available at github: https://github.com/hiljedoekes/PhageCom.

4) The parameters tested should be in the main text. Note that, other than the construction of the equations according to life history parameters, this is the only connection to "biology" that non-theoretical researchers will have to the manuscript. Table A1 should be in the main text. The table should list the conventional units of the parameters. For instance does burst size of 2 refer to 2 virions? This value would normally be in the 10's-100's for most phages. Some of these parameters don't seem reasonable at first glance. For instance, the adsorption rate is usually something on the order of 10-10 to 10-12. I suspect this could be because these are "scaled" values. If the scaled values are made for model convenience they should also be listed separately. This way the connection between the model and the biology will be more evident to readers.

Thank you for pointing out that including the parameter values in the main text will help non-theoretical researchers relate to the model. This is something we had not considered. We have moved the parameter table from the appendix to the Materials and methods.

Indeed, the parameters originally presented in Table A1 were scaled parameters (which are unit-less). The scaled burst size (or “effective burst size”, as we call it) is a composite of the actual burst size B and the probability b that if a phage adsorbs to a susceptible cell, that cell actually gets infected.

Consequently, it is indeed lower than commonly known values of the burst size. The scaled adsorption rate is a composite of the “plain” adsorption rate, the reproduction rate of bacteria (such that time is measured in bacterial generations, rather than hours), and the carrying capacity of bacteria. As a result, it is indeed much higher than the commonly known adsorption rates per minute.

We understand that presenting the scaled parameters only might cause confusion. We have therefore extended Table 1 to now include both the estimates for the non-scaled parameters that we took from the literature, and the corresponding values for the scaled parameters as they were used in the parameter sweeps. This is now also briefly explained in the Materials and methods.

5) The idea of having to go away from a steady state for this to work is very appealing. It is also likely much closer to biological reality in this case. Other than models that specifically address serial passaging, is the approach of disturbing ESS's regularly a new concept? Is there a generalizable framework for this?

Thank you for this comment, we agree that the need to consider the system away from steady state is one of the key points in our work. We are not familiar with any works that describe a general framework for determining ESS’s in disturbed systems. While we appreciate the idea of developing such a general framework, we are hesitant to overstate our case. We expect that there are many systems (both in models and in the natural world) in which disturbances away from steady state drive evolution, but the exact nature of these disturbances and the way in which they affect selection pressures might differ from system to system. It is hence hard to imagine a general framework that goes beyond merely stating that “disturbances should be included in ESS calculations”, and we are currently unsure what such a general framework should look like.

6) Some conclusions should be specified or clarified:

a) Figure 2 B-C: Some diversity seems to survive at steady state for T>5h. I believe this is what the authors refer to as "quasi-species" in the manuscript. It would be important to discuss and interpret this surviving diversity, and to explain whether it is going to survive forever or to slowly decay, and why.

This remaining diversity is indeed the viral quasi-species (Eigen, 1971). The quasispecies is shaped by the mutation-selection balance. While a single phage variant is the evolutionarily stable strategy (the most abundant variant, with φ = 0.04), phage variants with similar but slightly different φ-values constantly arise from mutations that occur during reproduction of the phages (i.e., during the production of new phage particles in the lytic cycle). Because these mutants only differ slightly from the ESS, they are likely also only slightly less fit. Hence, selection against them is quite weak, and these mutants will linger for some time, also reproducing themselves (which yields additional mutants, that differ more from the ESS and hence experience stronger selection). Some diversity in the quasi-species will thus be maintained forever, and the level of this diversity is determined by a balance between influx of mutants because of mutations arising during the reproduction of the ESS-phage variant (mutation) and decline of the density of these mutants because they are outcompeted by the ESS-phage (selection). Only if we would remove mutations from our model (i.e., set the mutation rate to zero) the diversity would decay until only the ESS-phage is left.

Since similar points were raised multiple times, we now realise that the quasi-species concept is not as widely known as we thought. We therefore included an explanation in our discussion of Figure 2B and C in the Results.

b) Figure 4: Please explain the difference between the analytical prediction and the simulation results at small values of bB in Figure 4B. In particular, does any of the simplifying assumptions result in this discrepancy?

The difference between the analytical prediction and the simulation at small values of bB indeed arises from violations of the simplifying assumptions we make to derive the ESS. The value of bB is the main determinant of the speed at which the epidemic unfolds. For small values of bB, the dynamics of the epidemic are slow. Then our simplifying assumption about the complete separation of time scales between the build-up of arbitrium during the epidemic and arbitrium uptake and degradation of by cells no longer holds. If the epidemic takes relatively long, the uptake and degradation of arbitrium by susceptible cells can no longer be ignored. As a consequence, the real arbitrium concentration at a certain time point (i.e., the moment phages “should” switch from a lytic to a lysogenic strategy) will be lower than the arbitrium concentration used in the analytical approximation. We should hence expect that the response thresholds found in the simulations are lower than the analytical predictions. This is indeed the case.

We included a short version of this explanation in the Results.

https://doi.org/10.7554/eLife.58410.sa2

Article and author information

Author details

  1. Hilje M Doekes

    1. Theoretical Biology, Department of Biology, Utrecht University, Utrecht, Netherlands
    2. Laboratory of Genetics, Department of Plant Sciences, Wageningen University, Wageningen, Netherlands
    Contribution
    Conceptualization, Software, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    For correspondence
    hiljedoekes@gmail.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6360-5176
  2. Glenn A Mulder

    Theoretical Biology, Department of Biology, Utrecht University, Utrecht, Netherlands
    Contribution
    Software, Investigation, Visualization, Writing - review and editing
    Competing interests
    No competing interests declared
  3. Rutger Hermsen

    Theoretical Biology, Department of Biology, Utrecht University, Utrecht, Netherlands
    Contribution
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4633-4877

Funding

Human Frontier Science Program (RGY0072/2015)

  • Hilje M Doekes
  • Rutger Hermsen

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

Acknowledgements

We thank Rob J de Boer for valuable discussions and comments on the manuscript. This work was supported by the Human Frontier Science Program, grant nr. RGY0072/2015 (http://www.hfsp.org/funding/research-grants).

Senior Editor

  1. Aleksandra M Walczak, École Normale Supérieure, France

Reviewing Editor

  1. Samuel L Díaz-Muñoz, University of California, Davis, United States

Reviewer

  1. Samuel L Díaz-Muñoz, University of California, Davis, United States

Publication history

  1. Received: April 29, 2020
  2. Accepted: January 15, 2021
  3. Accepted Manuscript published: January 18, 2021 (version 1)
  4. Version of Record published: March 5, 2021 (version 2)

Copyright

© 2021, Doekes 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

  • 1,432
    Page views
  • 187
    Downloads
  • 1
    Citations

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

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

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

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

  1. Further reading

Further reading

    1. Cancer Biology
    2. Computational and Systems Biology
    Kevin Hu et al.
    Research Article Updated

    In cancer, telomere maintenance is critical for the development of replicative immortality. Using genome sequences from the Cancer Cell Line Encyclopedia and Genomics of Drug Sensitivity in Cancer Project, we calculated telomere content across 1299 cancer cell lines. We find that telomerase reverse transcriptase (TERT) expression correlates with telomere content in lung, central nervous system, and leukemia cell lines. Using CRISPR/Cas9 screening data, we show that lower telomeric content is associated with dependency of CST telomere maintenance genes. Increased dependencies of shelterin members are associated with wild-type TP53 status. Investigating the epigenetic regulation of TERT, we find widespread allele-specific expression in promoter-wildtype contexts. TERT promoter-mutant cell lines exhibit hypomethylation at PRC2-repressed regions, suggesting a cooperative global epigenetic state in the reactivation of telomerase. By incorporating telomere content with genomic features across comprehensively characterized cell lines, we provide further insights into the role of telomere regulation in cancer immortality.

    1. Biochemistry and Chemical Biology
    2. Computational and Systems Biology
    Holly Fowle et al.
    Research Article

    Protein phosphorylation is a reversible post-translation modification essential in cell signaling. This study addresses a long-standing question as to how the most abundant serine/threonine Protein Phosphatase 2 (PP2A) holoenzyme, PP2A/B55α, specifically recognizes substrates and presents them to the enzyme active site. Here, we show how the PP2A regulatory subunit B55α recruits p107, a pRB-related tumor suppressor and B55α substrate. Using molecular and cellular approaches, we identified a conserved region 1 (R1, residues 615-626) encompassing the strongest p107 binding site. This enabled us to identify an 'HxRVxxV619-625' short linear motif (SLiM) in p107 as necessary for B55α binding and dephosphorylation of the proximal pSer-615 in vitro and in cells. Numerous B55α/PP2A substrates, including TAU, contain a related SLiM C-terminal from a proximal phosphosite, 'p[ST]-P-x(4,10)-[RK]-V-x-x-[VI]-R'. Mutation of conserved SLiM residues in TAU dramatically inhibits dephosphorylation by PP2A/B55α, validating its generality. A data-guided computational model details the interaction of residues from the conserved p107 SLiM, the B55α groove, and phosphosite presentation. Altogether these data provide key insights into PP2A/B55α mechanisms of substrate recruitment and active site engagement, and also facilitate identification and validation of new substrates, a key step towards understanding PP2A/B55α role in multiple cellular processes.