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).
https://doi.org/10.1371/journal.pcbi.1013340.g001
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).
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:
(1)
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.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 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 of 10μm, and volume Vol of 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 that monitor how many free phages of a specific type are present at any point in time and space. In each simulation time-step, , ‘bacterial growth’, ‘adsorption’, ‘infection’, ‘lysis’, ‘decay’ and ‘diffusion’ substeps occur (Fig 2).
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.
https://doi.org/10.1371/journal.pcbi.1013340.g002
Adsorption (Fig 2 step 1): the number of adsorption events for each phage strain is drawn from a Poisson distribution with mean 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 where and 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 .
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 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 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).
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, , used is Gaussian, and therefore symmetric. Asymmetric lysis time distributions are considered in Figs B and E in S1 Text.
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.
https://doi.org/10.1371/journal.pcbi.1013340.g003
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 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).
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 .
https://doi.org/10.1371/journal.pcbi.1013340.g004
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
(4)
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 , 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 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 : 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 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 ( ). 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 ( ). 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).
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 . 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.
https://doi.org/10.1371/journal.pcbi.1013340.g005
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.
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 . 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.
https://doi.org/10.1371/journal.pcbi.1013340.g006
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 m/min) than in the normal model (max approx. 11 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.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].
Live and let die: Lysis time variability and resource limitation shape lytic bacteriophage fitness
