Live and let die: Lysis time variability and resource limitation shape lytic bacteriophage fitness

TL;DR

This study shows that variability in bacteriophage lysis times can give them a competitive edge, even if it slows their growth, highlighting the importance of considering such variability in phage research.

Contribution

The study identifies a novel effect called 'population resonance' where lysis time variability enhances phage fitness in competitive scenarios.

Findings

Phage analytic growth rate in isolation poorly predicts fitness in competition experiments.

Lysis time variability provides a fitness advantage in serial passage and plaque expansion.

Models ignoring variability may miss viable evolutionary strategies in phage competition.

Abstract

Bacteriophages (phages) play a critical role in controlling bacterial populations, both in nature and as potential therapeutic agents. Their ability to replicate, compete against each other, and eradicate target cell populations is usually understood through a number of ‘life history parameters’, traditionally measured by population-level assays, which implicitly average the parameter’s value across a large number of infection events. Recent experiments suggest that bacteriophage life history parameters are subject to considerable stochasticity, raising the question of whether experimental and modelling efforts that do not account for this variability may overlook important factors in phage’s behaviour, competitive fitness or therapeutic viability. Here, using agent-based simulations, we investigate the importance of stochasticity in lysis time and burst size of lytic bacteriophages in two common laboratory competition experiments: serial passage of well-mixed populations and plaque expansion across a bacterial lawn. We find that a phage’s analytic growth rate in isolation can be a poor predictor of its fitness advantage in simulated competition experiments. Specifically, when lysis times are tightly distributed, we identify a novel effect we name “population resonance”, through which a bacteriophage can display a significant fitness advantage over a competitor with a much greater growth rate in isolation. Our simulations also show that both serial passage and plaque expansion reward variability in lysis time more than expected, by increasing the phage resilience when resources are scarce.

Bacteriophages (viruses that infect bacteria) can be described by a set of attributes called ‘life history parameters’. Historically, these parameters could only be measured on-average over large populations. Recent experimental advancements, however, have enabled their quantification at the single-cell-single-virus level, and have revealed that these traits are subject to inherent variability, even for genetically identical viruses infecting genetically identical cells. Here, we explore whether the degree of variability in two particular life history parameters; the lysis time (the time taken for a bacteriophage to infect and kill its host bacterium) and burst size (the number of new viruses produced per infected bacterium) might be subject to natural selection, by quantifying their effects on a phage’s fitness in two simulated competition experiments: serial passage and plaque expansion. Surprisingly, we find that lysis time variability is advantageous in competition even when it reduces the phage’s growth rate in isolation. This finding opens the door to future computational and experimental work, as it demonstrates that mean values alone are not sufficient to describe or predict a bacteriophage efficacy, and that models that ignore this variability can overlook viable evolutionary strategies.

Figures

a. Schematic of bacteriophage lytic cycle, characterised by seven life history parameters. b. Assembly curve, showing the number of viable virions within an infected cell as a function of time: initially zero until the eclipse period elapses, then increasing linearly. c to e. Assembly curves, each showing the effect of a faster-than-average lysis cycle (purple) versus the average (red), assuming the cycle is faster because the eclipse period is shorter (c), the post-eclipse is shorter (d), or both are proportionally shorter (e).

Schematic description of the agent-based model’s implementation.Healthy cells and progressing growth clocks are pictured yellow. Infected cells, viable virions and progressing lysis clocks are red. Virions that have injected their DNA, or that have decayed are pictured grey, as well as bacterial growth clocks which have halted due to infection.

Dependence of bacteriophage growth rate on lysis time under optimal conditions.a. Bacteriophage effective growth rate k as a function of mean lysis time, with either 0 minutes or 2.5 minutes standard deviation in lysis time. Inset: example lysis time distributions with mean 17 minutes, std 0 (solid line) or 2.5 minutes (dashed line). b. Optimal lysis time under ‘linear with eclipse’ model as a function of lysis time std. c. Effective growth rate achieved at optimal lysis time, under ‘linear with eclipse’ model, as a function of lysis time std. d. to f. Effective growth rate as a function of both mean lysis time and lysis time std for ‘normal’, ‘’linear‘, and ‘linear with eclipse’ models respectively.

Results of serial passage competition experiment, competing a wild type against a mutant phage with different mean lysis times.i. Difference in numerically predicted growth rate, using Eq 3. ii. Mutant fixation probability in stochastic simulation. iii. Fitness difference in stochastic simulation. Note that in the first and third rows, the visual scale has been set to saturate at ±0.1.

Results of serial passage competition experiment, in which a wild type phage with lysis time mean 17 min, std 2.5 min competes against a mutant with varied lysis time mean (y-axis) and std (x-axis).i. Difference in numerically predicted growth rate, using Eq 3. ii. Mutant fixation probability in stochastic simulation. iii. Fitness difference in stochastic simulation. Note that in the first and third rows, the visual scale has been set to saturate at ±0.1. Axes are scaled such that each pixel step in the y direction corresponds to a 1 minute change in mean lysis time, while each pixel step in the x direction corresponds to an equivalent fractional change in lysis time std versus the wild type value: 2.5/17 minutes. The four-point star on all panels correspond to the point in parameter space at which both phages have exactly the same parameters.

Plaque expansion simulation.a. Plaque expansion experiment schematic. b. Results of plaque expansion competition experiment, in which a wild type phage with lysis time mean 17 min, std 2.5 min competes against a mutant with varied lysis time mean (y-axis) and std (x-axis). i. Difference in numerically predicted growth rate, using Eq 3. ii. Mutant fixation probability in stochastic simulation. iii. Fitness difference in stochastic simulation. iv. Infection front velocity. Note that in the first and third rows, the visual scale has been set to saturate at ±0.05. Axes are scaled such that each pixel step in the y direction corresponds to a 1 minute change in mean lysis time, while each pixel step in the x direction corresponds to an equivalent fractional change in lysis time std versus the wild type value: 2.5/17 minutes. The four-point star on all panels corresponds to the point in parameter space at which both phages have exactly the same parameters. The white polygons in rows ii and iv outline the parameter space in which the wild type reaches fixation.

Equations

Full text

1 Introduction

Bacteriophages (phages) play a critical role in controlling bacterial populations both in nature, in which it is estimated that bacteriophages are responsible for 20–40% of all cell lysis events globally [1], and in medicine, where they provide a potential solution to the increasing spread of antimicrobial-resistance among human pathogens [2,3]. In fact, bacteriophages have been successfully used as antimicrobials since Felix d’Herell’s experiments on rabbits in 1917 [4,5], and recent trials have shown encouraging results [6,7].

Bacteriophages can be described by a number of ‘life history parameters’, including adsorption rate (how readily the phage binds to target cells), lysis time (the time between adsorption to the cell and lysis), and burst size (the number of offspring released at cell death). The phenotypic assays typically used to measure these parameters were first developed several decades ago, and rely on population-level measurements, which inevitably average a parameter’s value across a large number of infection events [8–10]. These assays have very recently seen some innovation to accommodate cell-to-cell variability [11–13], but in experiments focussed on phage evolution, competition and optimality, measuring the population average life history parameters remains the popular approach [14–19].

In parallel to this trend, most computational and theoretical work on bacteriophages has traditionally employed ordinary or partial differential equation methods, modelling the progression from infection to lysis either as an explicit and deterministic delay, [20–23], as an exponential process [24,25], or as a chain of intermediate steps connected by exponential processes [26]. These approaches assume that the distribution of lysis times will conform to a delta function, an exponential distribution, or an Erlang distribution, respectively. Some such work considers plasticity in viral parameters in response to the availability of nutrients to the bacterial host [27,28]. Accounting for such extrinsic variation allows for a deeper understanding of phage population and evolutionary dynamics, but the intrinsic variation that exists even between clonal cells in identical conditions [29] (and therefore between infection outcome) remains unexplored.

Agent-based models (ABMs), in which bacteria and/or phages are treated as discrete objects with stochastic behaviour, are particularly well suited to investigate the role of noise in population dynamics. Specifically for phage populations, they have been successfully used to study the evolutionary consequences of a phage permitting or preventing superinfection [30], the physical volume occupied by bacteria in a growing biofilm [31,32], the decay in heterozygosity at the front of an expanding plaque [33], and the possibility of a phage evolving towards symbiosis [34]. By construction, ABMs introduce a degree of variability, since the progression from each time step to the next is probabilistic, but few [30], to our knowledge, have incorporated explicit heterogeneity in lysis time or burst size.

In stark contrast to the general trend in the field summarized above, very recent experiments suggest that bacteriophages are subject to considerable heterogeneity in adsorption rate [16], lysis time [16,35,36] and burst size [12,16,37], the latter varying by more than an order of magnitude even in genetically identical cells in homogeneous conditions [37]. Since phages are totally dependent on the molecular machinery of their host, this variability likely arises from the considerable heterogeneity of the target bacteria’s internal metabolic state, via the availability of ribosomes [38–40], surface receptors [41,42] and other critical molecules. Adsorption rate has also been shown to be influenced by environmental conditions such as pH and temperature directly [43]. Experimental and modelling efforts that do not account for such variability may, therefore, overlook critical factors that influence phage’s behaviour, competitive fitness or therapeutic efficacy [44].

While the ultimate source of variability in life history parameters likely arises from the host cells, experimental evidence suggests that a bacteriophage can modulate this variability, as different genetic backgrounds exhibit it at different levels [35,36,45]. We may therefore think of a “noisy phage” as a phage that possesses some genetic mechanism to amplify the biological noise of its host, and vice versa, a “non-noisy phage” as one suppressing it. One genetic mechanism which has been shown, for instance, to influence lysis time variability is stochastic expression of holins, proteins which create gaps in the host bacteria’s surface membrane [35,36]. As long as these mechanisms are encoded in the phage genome, we can think of them as inheritable traits subject to selection.

In bacterial populations, variability in doubling time has been shown to result in a population growing faster than its mean doubling time would predict [46,47], but whether this result holds in phage populations is not obvious. A key distinction is that, while a bacterium always produces two offspring per division, bacteriophage burst sizes are much larger, can themselves be stochastic [12,37] and, crucially, may strongly depend on lysis time [19,23]. Here, we elucidate the importance of stochasticity in lysis time and burst size in lytic bacteriophage replication by using agent-based simulations that explicitly model lysis time and burst size distributions and their potential inter-dependencies in two common laboratory competition experiments: serial passage of well-mixed populations and plaque expansion across a bacterial lawn.

We find that a phage’s analytic growth rate in isolation is a poor predictor of the outcome of competition experiments and fails to capture two key features associated with lysis time variability. First, both the serial passage and plaque expansion protocols reward variability in lysis time more than expected: a “noisier” phage is able to outcompete a faster growing, less noisy phage. Second, when lysis times are tightly distributed, we observe an effect we name “population resonance”, in which one phage displays a significant fitness advantage over a competitor with a much greater growth rate in isolation. Both of these phenomena arise from the resource-limited nature of the phage environment and reveal distinct strategies that a phage can adopt to outlive its viral competitors.

2 Methods

2.1 Lytic bacteriophages can be described by seven life history parameters

The lytic cycle (Fig 1A) begins with a phage particle (virion) adsorbing to the surface of a susceptible cell. The virion then injects its genetic material (DNA or RNA, depending on the phage species), which hijacks the cell’s molecular machinery. The cell is forced to synthesise many copies of the phage’s proteins and genome, before finally producing lysis proteins, called holins and lysins, which cause the cell to lyse (burst), releasing the assembled virions [48,49]. Virions then seek new prey through passive diffusion, and may stochastically denature.

a. Schematic of bacteriophage lytic cycle, characterised by seven life history parameters. b. Assembly curve, showing the number of viable virions within an infected cell as a function of time: initially zero until the eclipse period elapses, then increasing linearly. c to e. Assembly curves, each showing the effect of a faster-than-average lysis cycle (purple) versus the average (red), assuming the cycle is faster because the eclipse period is shorter (c), the post-eclipse is shorter (d), or both are proportionally shorter (e).

To compare and predict the action of different bacteriophage strains, we employ a widely-used model [23,25,27,28,33,50] in which a bacteriophage is described by seven life history parameters (Table 1 and Fig 1A).

Table 1: Viral life history parameters, their units and definitions.

Importantly, these life history parameters do not vary independently, as physics and biology impose both interdependencies and absolute limits. For instance, to release more virions at lysis (greater burst size, β), a cell must undertake virion production for longer (greater lysis time, τ) or assemble virions faster (greater maturation rate, m). We model this dependency through Eq 1. This equation has previously been used to describe the assembly of T7 virions by prematurely and externally lysing host cells [19], as well as to relate the lysis times and burst sizes of a number of mutants derived from the same ancestral strain of bacteriophage λ [23]. In line with this body of work, we define the relationship between a bacteriophage’s burst size β and lysis time τ as:

$$[eqn]$$

According to this model, the number of viable virions within a cell increases linearly with time, after some initial eclipse period (Fig 1B). This assumption is consistent with experimental evidence [19,37] showing that viral components are produced almost linearly until lysis occurs, due to saturation of the host cell’s ribosomes. If lysis times were very long, we expect that this linear relationship might break down, with the assembly curve likely approaching a sigmoid, as the cell’s total available amino acids and nucleotides become the limiting factor (yield-limitation), rather than its ribosomes (rate-limitation). Since no experimental data is currently available in these regimes, here we will focus on the rate-limited case, though all the methods used here could be extended to a logistic model of virion assembly.

While in principle all three parameters relating burst size and lysis time, m, τ and ε, could evolve, it is typically assumed [23] that m and ε are set by some absolute biological constraints (for a specific bacterial genotype in specific environmental conditions): the production rate m is as fast as biologically possible and the eclipse time ε is as short as biologically possible. In other words, a larger value of m and smaller value of ε would always lead to a fitness advantage, so, if the phage could evolve to further increase m or decrease ε, it would already have done so. The possibility of cell-to-cell variation in m is discussed in Sect 2.2.

As a result, the only remaining parameter that can be subject to selection is τ: Eq 1 introduces a strategic decision between producing fewer offspring quickly or more offspring slowly. The consequences of these strategies on a bacteriophage’s population-level growth rate, as well as how variability in lysis time implies variability in burst size are discussed in the following section.

2.2 Variation in lysis time implies correlated variation in burst size

Assuming virions are produced according to Eq 1, there are several possibilities as to how stochasticity in lysis time (τ) may create correlated variation in burst size (β).

(i) If lysis time stochasticity arises purely from variation in the duration of the eclipse period and no variability in production time is observed, no correlated variation in burst size is expected (Fig 1C). Burst size may still vary from cell-to-cell owing to intrinsic stochasticity in maturation rate, but, crucially, burst size and lysis time will be uncorrelated at the single-cell-single-virion level. We label this possibility the ‘normal’ model. In particular, since in this model lysis time and burst size are uncorrelated, it also represents a null-model where no trade-off between lysis time and burst size is present.(ii) If lysis time stochasticity arises purely from variation in the post-eclipse period of active virion production, then given the lysis time of an individual virion, we can predict the burst size of that virion by applying Eq 1 directly (Fig 1D). We label this the ‘linear with eclipse’ model, or just ‘eclipse’ model.(iii) If lysis time stochasticity affects both the eclipse and post-eclipse proportionately, the relative variation in burst size will be equivalent to the relative variation in lysis time (e.g., 20% reduction in lysis time will lead to 20% reduction in burst size, Fig 1E). We label this the ‘linear’ model.

In the following, we consider each of these three models as a biological possibility, and determine in which cases variability in lysis time is advantageous. We focus on variability in lysis time specifically, as variability in burst size alone has little to no effect on a bacteriophage’s ability to propagate itself or clear a cell culture (Fig A in S1 Text).

2.3 Bacteriophage growth rate in isolation can be predicted by numerical integration

To provide a theoretical expectation for the effect of lysis time variability for the three different models discussed above, we consider a bacteriophage with lysis time distribution $[eqn]$ , and burst size $[eqn]$ , and predict its growth rate under the assumption that bacteria are abundant, and therefore: (a) adsorption is instantaneous, and (b) the probability of adsorbing to an already-infected cell is vanishingly small. Under these assumptions, the population of bacteriophages as a function of time ( $[eqn]$ ) is

$$[eqn]$$

where the effective growth rate k is given by

$$[eqn]$$

which can be computed numerically for any distribution $[eqn]$ . The extension of this result to a low-density bacterial culture in which adsorption is not fast relative to lysis time can be found in Fig K in S1 Text.

2.4 Agent-based simulation of two-phage competition experiments

The numerical calculation of the effective growth rate k derived above provides a measure of phage selective advantage under the assumption that resources (bacteria) are so abundant that phages do not meaningfully compete against each other. In real case scenarios, including laboratory experiments, phages always interact, if nothing else because of competition for resources. To assess the competitive fitness advantage of a phage population under different stochastic models, we implement two types of agent-based simulations, mimicking two distinct common laboratory protocols used in the study of bacteriophage evolution: serial passage (spatially unstructured) [14,51] and plaque expansion across a bacterial lawn (spatially structured) [52,53]. The former is characterized by a total well-mixed volume Vol of $[eqn]$ ml (carrying capacity 100,000 cells) initialized with 10,000 uninfected cells at random points in their cell cycle. The latter is modelled as a one-dimensional series of spatial ‘demes’, each corresponding to a length $[eqn]$ of 10μm, and volume Vol of $[eqn]$ ml (carrying capacity 100 cells), initialized with 100 uninfected cells at random points in their cell cycle. Note that we model a bacterial lawn, not a biofilm and as such all cells are considered to be metabolically active. While the volumes and cell counts used in both simulations are much lower than would be used in a typical laboratory experiment, our model and results are independent of absolute scale (Fig F in S1 Text), and the values have been chosen to reproduce experimentally observed adsorption curves in well-mixed cultures [8] and plaque expansion speeds on bacterial lawns [25,54]. We can therefore think of these simulations as representing a subsection of a larger experiment.

Each bacterial agent is characterized by two internal clocks, one representing its own cycle of growth and division and the other representing a count down to lysis, as well as an infection state: ‘healthy’, ‘infected (wild type)’ or ‘infected (mutant)’, and an initially empty list of adsorbed virions. Phage populations are modelled as counters $[eqn]$ that monitor how many free phages of a specific type are present at any point in time and space. In each simulation time-step, $[eqn]$ , ‘bacterial growth’, ‘adsorption’, ‘infection’, ‘lysis’, ‘decay’ and ‘diffusion’ substeps occur (Fig 2).

Schematic description of the agent-based model’s implementation.Healthy cells and progressing growth clocks are pictured yellow. Infected cells, viable virions and progressing lysis clocks are red. Virions that have injected their DNA, or that have decayed are pictured grey, as well as bacterial growth clocks which have halted due to infection.

Adsorption (Fig 2 step 1): the number of adsorption events for each phage strain is drawn from a Poisson distribution with mean $[eqn]$ Vol, where B is the local number of uninfected bacteria, I is the total local number of infected bacteria, and Vol is the local volume. A bacterial target, sampled randomly and with replacement, is then associated to each of these events. Adsorptions have a chance to fail proportional to the number of phages already adsorbed to the target cell, representing the finite number of receptor sites each cell possesses. The phage population counter is then updated depending on the number of successful adsorption events.

Infection (Fig 2 step 2): Uninfected bacteria that have at least one adsorbed phage move to an ‘infected (wild type)‘ state or ‘infected (mutant)’ state according to the nature of the adsorbed phages. In the case where a cell has adsorbed phages of both types, the type of infected state is chosen binomially using the proportion of adsorbed phages of each type. Once a cell becomes infected, the lysis time is drawn from a normal distribution $[eqn]$ where $[eqn]$ and $[eqn]$ correspond to the lysis time mean and std of the infecting phage, the lysis-clock starts, and the growth-clock is paused. We note that in this setup, there is no benefit to a bacteriophage that adsorbs to an already-infected cell at a later time step, since the first successful infection immediately and permanently determines the cell’s fate. We use the term “super-adsorption” to refer to a virion binding to an already-infected cell, as the secondary virion is removed from the pool but gains no benefit, since the status of the cell is unchanged. This contrasts the phenomenon of “super-infection” [30], in which the secondary virion is able to secure a share of cellular resources, and thus a share of the viral progeny.

Growth (Fig 2 step 3): uninfected bacteria advance their growth-clocks. If the growth clock completes its cycle, the cell replicates. A new growth cycle length is drawn for both mother and daughter cells from a normal distribution.

Lysis (Fig 2 step 4): infected cells advance their lysis-clocks. If the lysis clock completes its cycle, the cell is deleted, and a number of free phage of the appropriate type are added to their corresponding population counter. Depending on the stochasticity model, this number may be drawn from a distribution, or calculated from the total lysis time.

Decay (Fig 2 step 5): each phage strain’s population counter is reduced by a number drawn from a Poisson distribution with mean $[eqn]$ .

Diffusion (only for plaque expansion, Fig 2 step 6): the number of migrating phages from one deme to its neighbours is drawn from a Poisson distribution with mean $[eqn]$ and the corresponding phage population counters are appropriately updated.

Treadmill (only for plaque expansion, Fig 6A): to allow the infection front to propagate for longer without having to simulate a very large number of demes, we take a ‘sliding window’ approach. If the deme with the largest number of free virions (total across both pools) is more than 30 demes from the leftmost deme (300 $[eqn]$ m, or 50% of the total spatial extent X of the simulation), the leftmost deme is deleted, all other demes are shunted one position to the left, and a new deme is created at the right hand edge of the simulation.

Transfer (only for serial passage, Fig 4A): if the total number of bacteria falls below 100 (1% of its initial value), we perform a 1000X dilution, so that every remaining bacterium (infected or uninfected) and free phage has 1 in 1000 chance of being transferred to a fresh volume containing 10,000 new uninfected cells. The dilution factor of 1000 [14] and decision to transfer after lysis of the host culture (rather than at fixed time intervals) [51] are both consistent with contemporary evolution experiments.

Both protocols are initialised with 2 pools of 100 phage each, one labelled ‘wild type’, the other labelled ‘mutant’. In the ‘plaque expansion’ all phages are initially placed in the leftmost deme. The serial passage simulations run either until one phage pool outnumbers the other 70:30, at which point it is declared the winner, or for a maximum of 600 minutes. The plaque expansion simulations always run for 180 minutes. The parameter values used in both these simulations are as follows (Tables 2, 3, 4, and 5).

Table 2: Parameters used in serial passage simulations.

Table 3: Parameters used in plaque expansion simulations.

Table 4: Bacteria parameters used in all simulations.

Table 5: Bacteriophage parameters used in numerical analysis and simulations.

3 Results

3.1 Numerical calculations provide theoretical expectations for optimal lysis times

In Sect 2.3, we described how a bacteriophage’s effective growth rate k may be numerically predicted. In Fig 3A, we report this prediction as a function of mean lysis time for a standard deviation of either 0 or 2.5 minutes, respectively, for the three burst size models. Parameters were chosen to be similar to those of bacteriophage T7, and can be found in Table 5. An analysis of the sensitivity of our results to these choices can be found in Figs L and M in S1 Text. In the following: (i) the eclipse period ε was taken to be 9.5 minutes, (ii) under all three models, a lysis time τ of 17 minutes leads, by construction, to a burst size β of 150, and (iii) the lysis time distribution, $[eqn]$ , used is Gaussian, and therefore symmetric. Asymmetric lysis time distributions are considered in Figs B and E in S1 Text.

Dependence of bacteriophage growth rate on lysis time under optimal conditions.a. Bacteriophage effective growth rate k as a function of mean lysis time, with either 0 minutes or 2.5 minutes standard deviation in lysis time. Inset: example lysis time distributions with mean 17 minutes, std 0 (solid line) or 2.5 minutes (dashed line). b. Optimal lysis time under ‘linear with eclipse’ model as a function of lysis time std. c. Effective growth rate achieved at optimal lysis time, under ‘linear with eclipse’ model, as a function of lysis time std. d. to f. Effective growth rate as a function of both mean lysis time and lysis time std for ‘normal’, ‘’linear‘, and ‘linear with eclipse’ models respectively.

For the deterministic case (std = 0 min), the effective growth rate k decreases with increasing lysis time under the ‘normal’ and ‘linear’ models. By contrast, the ‘linear with eclipse’ model exhibits an optimum lysis time value due to the trade-off between long lysis times leading to very slow replication and short lysis time being unproductive because of the eclipse period. This preference for intermediate lysis times has been observed experimentally in serial passage experiments [14]. We emphasise that the fact that shorter lysis times lead to proportionately smaller bursts, as in the ‘linear’ model, is not sufficient to create an optimal lysis time under exponential growth conditions: the smaller burst size is always more than compensated for by the faster lysis cycle. The existence of this optimum depends entirely on a non-zero eclipse period ε.

When variability in lysis time is introduced, both the ‘normal’ and ‘linear’ models display an increase in growth rate compared to their deterministic counterparts. Mathematically, this feature can be explained by the convexity of the growth rate function: the mean of a symmetric sample about a given point will yield a greater value than the point itself. Importantly, the variation is symmetric in lysis time. Symmetric variation in lysis rate would cause a decrease in growth rate [57], since the function would no longer be convex. Conversely, the ‘linear with eclipse’ model, which features a concave function for short lysis times, is characterized by a reduction in growth rate with the optimal value shifting rightward (Fig 3B and 3C). Under this model, the noisier a phage’s lysis cycle is, the longer its lysis time should be and the lower the phage maximum growth rate is. This phenomenon holds across a large range of values, as shown in Fig 3D, 3E, and 3F.

Many other probability distributions can be used to investigate the role of noise. A biologically intuitive choice for lysis time distribution is the positively-skewed Erlang distribution, as it describes the stochasticity in waiting time from a series of n (shape) consecutive reactions, each occurring with rate λ (rate). An Erlang with shape $[eqn]$ is simply an exponential distribution. We find that the behaviour of the effective growth rate is qualitatively similar to that obtained using the Gaussian distribution (Fig B in S1 Text): ‘Normal’ and ‘Linear’ models favour short mean lysis time and large variance (low shape), while the ‘Linear with Eclipse’ model favours an intermediate optimal mean lysis time and small variance (high shape). A similar trend is also observed for the Weibull distribution (Fig C in S1 Text). To look specifically at the effect of asymmetry in the distribution, we have also considered a skew normal distribution (Figs D and E in S1 Text). Interestingly, we see that, for given mean lysis time, negative skew is beneficial in the ‘Normal’ and ‘linear’ models, as it provides additional weight to short lysis times, while positive skew is beneficial in the ‘eclipse’ model, as its short left tail limits the occurrence of very short unproductive lytic cycles. Given the qualitative similarities between the different distributions, in the following, we focus on the easily interpretable Gaussian distribution.

3.2 Serial passage competition experiments reveal winning strategies beyond growth rate maximisation

In Sect 2.4, we described a stochastic, agent-based model of the bacteriophage lytic cycle, which allows us to explicitly model the consequences of variability in a bacteriophage’s burst size and lysis time. We start by using this model to investigate how the phage fitness landscape is affected by (i) the different types of trade-off between lysis time and burst size encoded in our models (Fig 1C, 1D, and 1E), and (ii) the level of noise in lysis time. We consider two bacteriophage populations, labelled “wild type” and “mutant”, and assign them both the same burst size model (‘normal’, ‘linear’, or ‘linear with eclipse’) and the same standard deviation in lysis time (0.01 or 2.5 minutes). We then vary the mean lysis time of each between 10 and 30 minutes in one-minute increments, and compete the two phage populations against each other through a serial passage assay (Fig 4A, Sect 2.4, initial MOI 0.01 per phage species, 0.02 total).

Results of serial passage competition experiment, competing a wild type against a mutant phage with different mean lysis times.i. Difference in numerically predicted growth rate, using Eq 3. ii. Mutant fixation probability in stochastic simulation. iii. Fitness difference in stochastic simulation. Note that in the first and third rows, the visual scale has been set to saturate at ±0.1.

For each set of simulation parameters, we quantify competitive fitness using three alternative measures: (i) the numerically predicted difference in effective growth rate k (defined in Eq 3) between the two phages (Fig 4BI); (ii) the mutant fixation probability defined as the fraction of twenty independent runs in which the mutant fixes in the population (see Sect 2.4, Fig 4 BII); (iii) the rate of fixation of the mutant (positive values) or wild-type (negative values) strain (Fig 4BIII), defined as

$$[eqn]$$

$[eqn]$ simplifies to the numerical difference in growth rates k in the limiting case where resources are abundant and the two phage populations are growing exponentially (derivation in S1 Text Sect 2).

In every case red corresponds to the wild type displaying a fitness advantage and blue to the mutant displaying a fitness advantage. By construction, the plots are antisymmetric about $[eqn]$ , therefore in the discussion below we will focus on the bottom-right half of the diagram without loss of generality.

Starting with the leftmost column, corresponding to the ‘normal’ model with narrow lysis time distribution (0.01 minutes standard deviation), the top row (numerically predicted growth rate) replicates the results in Fig 3 showing that short lysis times lead to a fitness advantage. The middle row (mutant fixation probability) corroborates this numerical result. The bottom row, (mutant fixation rate), reveals an interesting finding. First, we see that while in the majority of parameter space the rate of fixation agrees in sign to the relative growth rate (top row), the rate of fixation increases much more slowly as we move away from the $[eqn]$ line. This result reflects how the competition for resources (unaccounted for in the relative growth rate) dramatically screens the fitness advantage that a fast lysing phage would have when grown in isolation (Figs J and N in S1 Text). A similar result was also observed in ref [30] and highlights how phage fitness estimates based on growth in isolation can be misleading and overestimate the fitness (dis)advantage of a population.

The second unexpected result is the appearance of a “stripe” corresponding to intermediate values (green) of fitness advantage within an otherwise highly advantageous region of parameter space (blue). Along this stripe, the mutant lyses cells faster than the wild type and reaches fixation as we would expect. However, it reaches fixation more slowly (in more serial passage cycles) than in the surrounding parameter space.

Close investigation of these cases reveals that they arise due to the limited number of cells available (Fig H in S1 Text). The faster-growing mutant experiences a substantial “dead time” after completing its final round of lysis, just before the transfer step takes place. During this period, the only cells remaining are those infected by the (slower-lysing) wild type. The mutant virions, therefore, have no way to further proliferate, and are only subject to decay or unproductive super-adsorption (“super-adsorption” rather than “super-infection”, as the secondary virion is removed from the pool but gains no benefit, as the status of the cell is unchanged). As a result, by the time the wild type finishes its final round of lysis, there are fewer mutant virions left to survive the dilution.

In line with this explanation, the lower edge of the stripe, highlighted with a white line, corresponds to $[eqn]$ : the longest possible dead time (and thus the least advantageous parameters for the mutant phage in this region) would occur when the end of the mutant’s final round of lysis happens just after the start of the wild type’s final round. When the mutant’s lysis time becomes half of that of the wild type, it will be able to begin a new round of lysis at the same time as the wild type’s final round, with minimal dead time (Fig H in S1 Text).

As additional support to our explanation, we repeat all simulations in Fig 4B but under “conservative” conditions, in which free virions do not decay and super-adsorption cannot occur (Fig J in S1 Text). As expected, under these conditions the stripe disappears, as the mutant phage suffers no penalty from spending “dead time” outside the cell. We refer to this phenomenon as “population resonance”, since it occurs only for specific ratios $[eqn]$ whose values depend on the number of lytic cycles between dilutions, which are themselves set by the initial multiplicity of infection (MOI) (Fig N in S1 Text).

The second column (‘normal’ model, 2.5 mins lysis time std) qualitatively agrees with the first. The only quantitative difference is that the stripe caused by population resonance appears blurred. Since here lysis times for both phages are more variable, lysing events are less synchronised, and the resonance effect is weaker.

The third and fourth columns of Fig 4B correspond to the ‘linear’ model ( $[eqn]$ ). When compared to the normal model, the principle difference is that the population resonance effect is stronger. The fixation probability plots (second row) show that, in the 0.01 min std case, population resonance is strong enough to flip which phage reaches fixation. In the 2.5 min case this is not true, as de-synchronisation in lysis cycles weakens the resonance effect, but slight deviations from 100% fixation are visible as a pale blue line in the otherwise uniform block of blue.

The reason why population resonance is stronger under the ‘linear’ model than under the ‘normal’ model is that in the ‘linear’ model longer lysis times confer larger bursts. Hence, while the faster-growing mutant is losing virions during the dead period, the slower-growing wild type benefits from a larger final burst. This longer lysis cycle and larger burst size would not be a winning strategy under exponential growth conditions, as shown in the theoretical expectation for the effective growth rate (top row). However, in the resource-limited conditions of a serial passage, forcing an opponent to lose virions to super-adsorption and decay while taking time to produce more offspring can be a winning strategy. This is in effect a weaponisation of a long lysis time.

Finally, the two rightmost columns on Fig 4B correspond to the ‘linear with eclipse’ model ( $[eqn]$ ). The topography of the fitness landscape is altered by the introduction of a critical point corresponding to the optimum lysis time, visible in the top row as the intersection of two (yellow) contours of zero fitness difference and corresponding to the maxima identified in Fig 3: at lower lysis time values in the low std case, and shifted to larger values in the higher std case. At the same time, the population resonance effect becomes even stronger than under the ‘linear’ model, as longer lysis times imply even greater benefits in burst size.

The fixation probability and the fixation rate plots for the 2.5 min std case (middle and bottom rows, rightmost column) shows that the population resonance stripes are now impactful enough to remain present even with greater lysis time variation, and hence less synchronised lysis cycles, but are “rounded off” by the noise in lysis time. Indeed, a very similar shape can be obtained by applying a Gaussian blur to the corresponding plots in the 0.01 min std case.

In summary, our results reveal the presence of population resonance as a novel effect that arises when two bacteriophage populations in competitions are characterized by specific relative lysis times. It is strongest when lysis time is tightly distributed, and stronger the more the trade-off between lysis time and burst size penalizes short lysis times. According to our simulations, population resonance allows one bacteriophage to outcompete another phage that would grow faster in isolation. One factor which makes population resonance easier to observe is our decision to carry out a transfer once 99% of the initial bacteria have been lysed (as in [51]), as opposed to after some fixed time interval (as in [19]). If a fixed interval were used, and this interval were longer than the time taken to lyse all cells in the culture, the only difference would be an extended period of virion decay between the lysis of the last cell and the transfer step. If both the wild type and mutant have the same decay rate (as in our case), the relative proportion of wild-type and mutant phage would not change, and the resonance is expected to hold in the same position. If, however, the transfer were carried out before the majority of cells had been lysed, additional resonance effects would occur, corresponding to fitness discontinuities as a phage’s lysis time aligns with specific fractions of the dilution interval, producing a more complex fitness landscape (Fig I in S1 Text).

3.3 Serial passage experiments suggest lysis time noise is beneficial even at the cost of reduced growth rate

Having shown that noise in lysis time meaningfully affects which of two competing bacteriophage strains reaches fixation, we now explore whether variability in lysis time can be under selective pressure and quantify its direct effects on competitive fitness. We repeat the serial passage competition experiment described in Sect 3.2, keeping the wild type phage’s mean lysis time and std constant (17 min and 2.5 min, respectively), while varying those of the mutant (Fig 5).

The first two columns, corresponding to the ‘normal’ and ‘linear’ models, are qualitatively similar to each other. In both cases, the theoretical growth rate difference (top row) suggests that a mutant with a shorter lysis time should grow faster than the wild type. The yellow contour of equal fitness has a slightly positive gradient, meaning that a phage should see an equal gain in fitness by either reducing its mean lysis time by 1 minute, or increasing its lysis time std by approximately 3 minutes.

The middle row (fixation probability) agrees with the overall claim that shorter lysis times are advantageous, but the slope of the equal-fitness contour is greater than in the top row, suggesting that, under the resource-limited conditions of the serial passage, each 1 minute decrease in lysis time mean is competitively equivalent to an increase of approximately 1 min in lysis time std. Since this increase in gradient similarly holds for both the ‘normal’ and the ‘linear’ models, it cannot be understood in terms of the trade-off between lysis time and burst size characteristic of the latter, but it is the manifestation of an intrinsic benefit of stochasticity under resource-limited conditions.

An explanation for this result is the following. In exponential growth, one bacteriophage strain will outcompete the other if its growth rate is consistently higher, therefore its lysis cycles must be fast on average. In the serial passage experiment, bacteria are limited and vastly outnumbered by the phage towards the end of each passage (MOI approaching 100). Here, a bacteriophage strain does not need to have short lysis cycles on average, rather it needs enough short lysis cycles to secure a larger share of the available bacteria, which can be achieved with large variation around the mean.

The rate of fixation (bottom row) provides a finer quantification of these results, as we see that the mutant fitness disadvantage at high mean lysis times predicted by growth in isolation can be overcome by increasing lysis time std. For a given std, the fitness difference as a function of mean lysis time displays a non-monotic behaviour. This is a result of the population resonance effect described in the previous section. Importantly, while increasing noise decreases the relative advantage emerging from the resonance effect, it also broadens the region of parameter space in which it is relevant (Fig O in S1 Text).

The rightmost column of Fig 5 corresponds to the ‘linear with eclipse‘ model. Here, the numerical growth rate prediction (top row) points to a small optimal region of parameter space (short lysis time, low std, lime green) below and to the left of the equivalence point, in which the mutant should grow faster than the wild type. This reflects the finding in Sect 2.3, and the fact that the wild type phage is already close to the optimum. Interestingly, the mutant fixation probability plot (second row) shows a different picture, in which two distinct regions of parameter space result in mutant fixation. The region towards the top of the plot is the result of population resonance (when super-adsorption and decay are disabled this region disappears, Fig P in S1 Text).

Results of serial passage competition experiment, in which a wild type phage with lysis time mean 17 min, std 2.5 min competes against a mutant with varied lysis time mean (y-axis) and std (x-axis).i. Difference in numerically predicted growth rate, using Eq 3. ii. Mutant fixation probability in stochastic simulation. iii. Fitness difference in stochastic simulation. Note that in the first and third rows, the visual scale has been set to saturate at ±0.1. Axes are scaled such that each pixel step in the y direction corresponds to a 1 minute change in mean lysis time, while each pixel step in the x direction corresponds to an equivalent fractional change in lysis time std versus the wild type value: 2.5/17 minutes. The four-point star on all panels correspond to the point in parameter space at which both phages have exactly the same parameters.

The high-fitness region at the bottom of the plot, which in the other models extends to the x-axis, is here instead bounded from below by the presence of the eclipse period: if the mean lysis time of the mutant is too close to the eclipse period, the mutant will not be competitive. This result agrees with numerical prediction (top row). The fitness difference in simulations, however, disagrees with numerical predictions in that the majority of this high fitness region is to the right of the equivalence point (larger stds), rather than below and to the left. This is because, as mentioned previously, when resources are sufficiently limited that only a small number of lysis cycles are possible, the optimal strategy is not to maximise steady-state growth rate, but rather to carry out fast early lysis cycles, which secure a majority share of the available resources, followed by slower but more productive later lysis cycles. An intermediate mean lysis time and larger lysis time std make this possible. The lysis time std cannot be made too great, or the phage risks lysing before or shortly after the eclipse period, producing nothing. Hence, the lower edge of this high-fitness region is slightly upward sloped, favoring longer mean lysis times when the std is large.

In summary, we find that serial passage rewards lysis time variability more than numerical growth rate predictions would suggest (Fig O in S1 Text), including some cases in which numerical growth rate wrongly predicts that variability would be detrimental. We reason that this is due to fast lysis cycles being disproportionately advantageous under resource-limited conditions as they ensure a greater share of the total, limited, pool of resources (bacteria). In the ‘linear’ and ‘linear with eclipse’ models, larger lysis time std allows for both shorter lysis cycles that produce virions quickly to claim the available resources, and longer, more productive later cycles that generate more virions. It should be stressed that this winning strategy is achieved only probabilistically: an individual phage is equally likely to cause a faster-than-average first lysis cycle as to cause a slower-than-average cycle. If the phage could sense the culture MOI and ‘decide’ to lyse rapidly when MOI is low and slowly when MOI is high, such as through a lysis-inhibition mechanism [58], this would be an even more competitive strategy.

3.4 Plaque expansion assay does not exhibit population resonance and results in weaker selective pressure on lysis time

We now move to a second standard phage growth assay, a plaque expansion across a bacterial lawn, which we simulate as a one dimensional series of well-mixed demes (Fig 6), as detailed in the Methods section. Free virions, but not bacteria, are able to diffuse between adjacent demes (Fig 6A). The initial MOI (per deme) is 1 per phage species, 2 total, in the leftmost deme, and 0 elsewhere, representing the typical local densities of viral particles in an initial inoculation of a clearing.

Plaque expansion simulation.a. Plaque expansion experiment schematic. b. Results of plaque expansion competition experiment, in which a wild type phage with lysis time mean 17 min, std 2.5 min competes against a mutant with varied lysis time mean (y-axis) and std (x-axis). i. Difference in numerically predicted growth rate, using Eq 3. ii. Mutant fixation probability in stochastic simulation. iii. Fitness difference in stochastic simulation. iv. Infection front velocity. Note that in the first and third rows, the visual scale has been set to saturate at ±0.05. Axes are scaled such that each pixel step in the y direction corresponds to a 1 minute change in mean lysis time, while each pixel step in the x direction corresponds to an equivalent fractional change in lysis time std versus the wild type value: 2.5/17 minutes. The four-point star on all panels corresponds to the point in parameter space at which both phages have exactly the same parameters. The white polygons in rows ii and iv outline the parameter space in which the wild type reaches fixation.

The data in the top row is identical to that in Fig 5, with a different colour range to make it directly comparable with the data from the simulations of plaque expansions. For the ‘normal’ and ‘linear’ models, the results share some common features with the serial passage case. Particularly, the line of equal fitness has a steeper gradient in the stochastic simulation plot (second row) than the predicted growth rate plot (top row), suggesting that lysis time variation is still rewarded more strongly than growth rate alone would predict. This is not surprising, as the typical MOI in plaque expansion is very high, making resources even more limited.

The fitness difference plots (third row) are markedly different to those produced by serial passage. First, the population resonance effect is no longer present, since virions no longer experience a “dead period” waiting for new bacteria to be added and lysis cycles quickly become desynchronised because of diffusion. Secondly, fitness (dis)advantages are overall smaller in magnitude here than in the serial passage case. In particular, a mutant with a lysis time mean larger than the wild type’s is able to persist for much longer and, vice-versa, a fast lysing mutant will take more time to fix.

The final row shows the infection wave’s velocity, which is determined by the winning phage. A mutant which produces a much faster infection front will quickly overtake and outcompete the wild type phage. The white polygon corresponds to the region of parameter space in which the wild type phage reaches fixation, and therefore we expect the wave velocity to be constant. Mutants with shorter mean lysis times and greater lysis time std produce faster infection waves, but infection waves are significantly slower in the linear model (max approx. 6 $[eqn]$ m/min) than in the normal model (max approx. 11 $[eqn]$ m/min), as a result of the trade-off between lysis time and burst size.

The third column corresponds to the ‘linear with eclipse’ model. As in the serial passage experiment, the optimal region of parameter space in which the mutant is able to reach fixation (second row, green ellipse) is not located where the numerical growth rate would suggest (top row). Its shift towards high std values strongly supports the idea that lysis time noise is beneficial in resource-limited environments.

The third row (fitness difference) and fourth row (infection wave velocity) show a general loss of contrast across the majority of the parameter space. This suggests that when a strong trade-off between lysis time and burst size is present (second and third columns), the plaque expansion puts little (second column, ‘linear’ model) to almost no (third column, ‘linear with eclipse’ model) selective pressure on lysis time, regardless of its distribution, as any benefit coming from lysing faster is compensated by the lower burst size.

In summary, a spatially structured environment removes any population resonance effect, and masks the selective pressure on lysis time compared to a serial passage setup. In line with this, we also find that, while noise in lysis time still confers a detectable fitness advantage, its overall impact is less significant, especially when a lysis time-burst size trade-off is present.

4 Discussion

4.1 Variability in lysis time as an evolvable trait with significant impact on fitness

Clonal bacteria growing under identical conditions display significant heterogeneity in gene expression due to intrinsic biochemical noise [29]. It is therefore expected that bacteriophages, which heavily rely on their host’s molecular machinery, are subject to similar unavoidable stochasticity in their replication. The extent and pervasiveness of such noise in a phage’s life history parameters cannot be measured by traditional phenotypic assays based on bulk population measurements, and as such, they have remained inaccessible until recently [12,36,37]. As a result, the vast majority of experimental and theoretical work aimed at quantifying phage ecology and evolution has focused on fitness effects linked to average values. Recent advances in single-cell measurements of phage-bacteria interactions [12,37] and more careful analyses of population assays [11] are finally enabling an experimental quantification of variability in life history parameters, inspiring parallel theoretical work.

Far from “noise” being simply an effect to which bacteriophages are subject, the degree of variability in a phage’s life history parameters has recently been shown to be, at least partially, under the phage genetic control [35,36,45], raising questions as to when and why such noise might be favoured by natural selection, or, alternatively, selected against. This parallels a broader trend in evolutionary biology, which sees stochasticity as not just influencing evolution, but, in some cases, driving it [59,60].

In this work, using extensive agent-based simulations, we have investigated how phenotypic stochasticity in lysis time and burst size can confer fitness (dis)advantages to a phage population in two commonly used experimental assays: serial passage [14] and plaque expansion [61]. Overall, we observe two surprising phenomena. First, variability in lysis time is generally an advantageous trait, even when it reduces growth rate in isolation. Second, when lysis cycles remain tightly synchronised, a novel effect we name “population resonance” may occur, in which a slower growing phage effectively weaponises a longer lysis cycle to create a sink for competing phages.

4.2 Lysis time variability is more advantageous than expected and its effects depend on the lysis time-burst size trade-off

One of the fundamental differences between noise in division time in bacteria and noise in lysis time in phages is the trade-off between the length of the latent period τ and the burst size β (unlike for bacteria, which will always replicate into two daughter cells upon division). For bacteria, it is known that variability in doubling time is advantageous, as long as it is symmetric [46]. The same is true for single-celled eukaryotes such as budding yeast [62]. By contrast, our results show that, for phages, the details of the trade-off between lysis time and burst size have strong influence on the fitness landscape and on how noise modulates it. Surprisingly, a linear correlation between burst size and lysis time results in a fitness map that is very similar to the one generated in the absence of any trade-off in both serial passage and plaque expansion scenarios. In both situations, we find a strong selection towards fast lysis cycles, whether this occurs by decreasing the mean lysis time, or by increasing the variability in lysis time so that more infection cycles are completed rapidly. For a non-zero optimal lysis time to emerge, a non-productive eclipse period has to be introduced in the trade-off. In this case, theory predicts that an optimal lysis time with minimal variation will produce the fastest growth rate. Interestingly, we find that, in both the serial passage and a plaque expansion competitions, variability in lysis time remains advantageous, even where it reduces population growth rate.

We find that the primary selector for such phenotypic variability is resource (host) limitation. Most laboratory settings, and likely wild environments, are characterized by periods of time or regions of space with very high MOI, where many phage virions compete for few host cells. In these conditions, noise in lysis time appears to confer increased resilience to a phage population. By producing virions throughout time and not in synchronized large bursts, variable lysis times can establish a more steady flow of virions, so that susceptible cells can be infected whenever and wherever they become available. This effect becomes crucial when free virions are subject to rapid loss, e.g., because of super-adsorption. Importantly, the desynchronisation achieved by lysis time variability is not equivalent to that achieved by a lower adsorption rate: a low adsorption rate inevitably reduces growth rate and competitive advantage because it increases the time spent by the phage searching for a host (see predicted growth rates in Fig 4 vs Fig K in S1 Text). By contrast, a variable lysis time with high adsorption rate allows any released phage to immediately infect a host a soon as it becomes available. Reduced adsorption rate may still be beneficial if prolonged periods of very high MOI are common, since it reduces a bacteriophage’s chance to be wasted through super-adsorption, but it would come at the cost of slowing down infection when susceptible cells do become available.

A loose analogy can be drawn between variability in lysis time and the stochastic emergence of bacterial ‘persisters’, the fraction of bacterial cells which enter a state of near-zero growth, even when resources are abundant [63]. Both strategies reduce the population’s growth rate while conditions remain ideal, but make it more resilient to massive loss under stressful conditions: extended periods of high MOI, or the application of antibiotics respectively.

4.3 Population resonance increases the number of viable lytic strategies

Our second key observation is the novel phenomenon we name ‘population resonance’, in which a phage which grows slowly in isolation is able to outcompete a faster-growing phage in serial passage experiments. In order for population resonance to occur, three conditions are necessary. First, lysis cycles must be comparatively synchronised. This means variation in lysis time must be small, and that adsorption must occur quickly, which implies that the cell culture must be relatively dense. Second, a trade-off between lysis time and burst size must be present, so that longer lysis times lead to larger burst sizes. Third, free virions must be subject to some decay, for instance via super-adsorption to previously infected cells. Since resources (susceptible bacteria) are limited, each phage will undergo a finite number of lysis cycles before the cell population is exhausted. If the slower-lysing phage’s lysis time value is just less than a multiple of the faster-lysing phage’s lysis time, for example 19 and 10 minutes respectively, it will be able to adsorb to the last uninfected cells present, forcing the final burst of fast-lysing phage to undergo a long “dead period” in which the free virions are subject to decay. By the time the slower-lysing phage completes its final lysis cycle, its offspring will outnumber what is left of the fast-lysing phages.

Combining our understanding of population resonance with the competitive advantage offered by lysis time noise provides an interpretation for the complex fitness landscape created by the serial passage experiment, particularly in the case of the ‘linear with eclipse’ model. We find that this assay favours two opposite strategies: “patient” and “rapacious”. The “patient” strategy uses a relatively slow mean lysis time with little variability, so that the phage can exploit the population resonance effect and effectively “outlive” the competing faster-lysing phage. The “rapacious” strategy favours shorter lysis time with high variability, so that the phage quickly removes potential host cells for the slower-lysing competing phage and relies on the occasional longer cycles for productive infections. Both the “rapacious” and the “patient” strategies are able to outcompete a “prudent” strategy that hinges on maximising growth rate in isolation (intermediate mean lysis time, very small std) because of the competition for limited resources (cells) between two phage populations. The efficacy of a similarly rapacious strategy was observed in [64], in which a fast-migrating, fast-lysing phage was able to out-compete a slower-lysing prudent phage, which grew more quickly in isolation. Conversely, lysis-inhibition [58] and lysogeny [65,66] which are both extant strategies in wild type phages, may be thought of as extensions of the “patient” strategy: when MOI is high, bursts are delayed as released virions would be likely lost to unproductive-super adsorption.

It would be interesting to test experimentally which strategy would actually preferentially emerge in an evolutionary experiment, since the serial passage is a well established laboratory protocol. While the outcome will be sensitive to the availability of mutations and potential epistasis, we expect that the rapacious strategy (short lysis, large variability) will typically be more robust, since the presence and location of population resonance, which is leveraged by the patient strategy, strongly depends on the details of the competing phage, the number of cells at the beginning of each serial passage, and the initial MOI (see Figs J, N, and P in S1 Text).

Turning to the plaque expansion assay, unexpectedly, we find that introducing an eclipse period in the lysis time-burst size trade-off almost completely flattens the fitness landscape, as long as lysis time is larger than the eclipse period. This result implies that there should be almost no selective pressure on lysis time or its variability in an evolutionary experiment performed on a bacterial lawn, possibly explaining why several studies performing such experiments predominantly find mutations related to adsorption [64,67]. Interestingly, our finding also suggests an indirect way to identify the presence of a strong trade-off between burst size and lysis time. A plaque assay competition dominated by fast-lysing phages would indicate the absence of a strong trade-off. In contrast, a plaque competition dominated by slow, highly productive phages could indicate the presence of a strong trade-off between lysis time and burst size, and, likely, the presence of a strongly conserved eclipse period. It is important to highlight that the plaque expansion assay is distinct from the infection of a bacterial biofilm. In plaque expansion, phages propagate across a metabolically homogenous bacterial lawn. By contrast, in a bacterial biofilm we would expect significant heterogeneity in the metabolic states of the bacteria [68,69]. Extending the model presented here to a bacterial biofilm would need to account for heterogeneity in nutrient availability, and therefore metabolic rate, as well as the presence of distinct bacterial phenotypes [68,69].

4.4 Population resonance has implications in both laboratory and wild environments

Population resonance, leveraged by the “patient” viral strategy, depends on synchronised infection events. In the serial passage protocol, this is achieved by the synchronous initial inoculation of phages into bacterial culture, followed by the synchronous inoculation into a fresh bacterial culture. While this synchrony is undoubtedly artificially amplified in some laboratory settings, population resonance may have broader applicability to other natural or artificial systems which are periodically flushed and in which bacteria and phages coexist, such as tidal aquatic environments [70], or the or human urinary [71] and gastrointestinal tract [72,73], if not as the primary driver of selection, then as an important ‘second-order’ consideration. For example, if one wishes to eradicate a pathogenic bacterial strain from the human GI or urinary tract, a cocktail of “rapacious” virulent phages and “patient” phages, potentially either lysogenic or possessing a lysis-inhibition-mechanism, might lead to better outcome. Reciprocally, if population resonance is predominantly confined to a laboratory setting, our results warn against over-interpretation of the emergence of optimal lysis times in experimental evolution, which might simply be the result of the selection protocol employed.

4.5 Summary

Predicting how an organism will fare in competition, based on its behaviour in an ideal environment with abundant resources and no meaningful competition, is challenging [74]. Resource investments ranging from toxin production in bacteria [75,76] to predation-defences in plants [77,78] reduce growth rate in monoculture, but become vital once competitors or predators are introduced. This phenomenon applies equally to bacteriophage. While phages lack the toxins and trichomes necessary to directly remove their competitors, we have demonstrated here that fast, noisy, or in some cases very slow lysis cycles can be weaponised to allow one phage to outcompete another, at the cost of reducing its own growth rate in isolation.

Can variability in lysis time be an adapted trait? When no or minimal trade-off between lysis time and burst size exists, our results show that selection will simply push the phage to lyse as quickly as possible. In this case, noisy phages could potentially emerge only if there is a lower limit to how short a deterministic lysis time can be or if it is evolutionary easier, because of availability of mutations or presence of epistatic effects, to increase variance in lysis time rather than decreasing its mean. However, in the presence of an eclipse period, an optimal lysis time value exists. Once the phage reaches this optimum, we expect that the selective pressure will concentrate on optimizing variability in lysis time. Experimentally testing these hypotheses will not only clarify the evolutionary role of noise in phage ecology and evolution, but also reveal which genetic changes can modulate variability in lysis time. Thus far, identified mutations which contribute to lysis time variability seem to hit the genes responsible for the expression of holin proteins, both in λ phage [36] and $[eqn]$ [45]. As holins are usually among the final genes expressed in the lytic cycle, partially responsible for lysis of the host [79], we expect that noise should affect only the post eclipse period, since virion production is already underway by the time holins are transcribed. Single-cell microfluidic experiments interestingly seem to suggest that the period between the beginning of capsid protein synthesis and lysis is under comparatively tight control, and that the period between adsorption and capsid expression is much more variable in wild-type phage T7 [37]. Taken together, these data suggests that variability in the lytic cycle of a bacteriophage is likely a complex phenomenon with no single source. Further experiments and genetic analysis of ‘noisy’ mutants are needed to deepen our understanding of the intricate interplay between an infecting phage and its host.

Supporting information

S1 TextAll supporting information is supplied in the material.(PDF)