Our system consists of two variants of the same phage (with population densities P₁ and P₂) that are identical in all respects except their propensity to go lysogenic. Both compete for a single bacterial strain in our model, but we separately track the density of uninfected bacteria (B₀), the density of bacteria infected by one or the other phage (B₁, B₂) and the density of lysogenized bacteria of each type (L₁, L₂). As we are modeling a well-mixed system, we use ordinary differential equations to represent the dynamics of the bacterial and phage populations:

where i can be 1 or 2, and B_tot and P_tot are the total bacteria and phage densities, respectively. This is a straightforward extension of the model Maynard et al. (2010) use to describe the dynamics of λ phage infecting a population of E. coli. The same equations have been used many times previously to model phage-bacteria interaction (Bohannan and Lenski, 2000; Beretta and Kuang, 2001; Weitz and Dushoff, 2007; Maynard et al., 2010; Sneppen et al., 2015). Parameters relating to bacteria are: γ, the bacterial growth rate, and K, the carrying capacity, i.e., the maximal total bacterial density the system can sustain. The other parameters are properties of the phages: the “burst size,” β, is the number of phage produced from the lytic infection of one bacterium; the “propensity of lysogeny,” f, is the fraction of infected bacteria that go lysogenic; η is the rate constant characterizing infection of bacteria by the phages.

The model assumes that: (i) “super-infection” of lysogens by either phage variant does not change the state of the cell (i.e., lysogens are immune to a subsequent infection), and (ii) “cross-infection” of infected bacteria by the other phage variant also does not change the state of the cell (i.e., the ongoing infection continues, and if the outcome is lytic, only the first type of phage is produced, while if the outcome is lysogeny, only the first type of phage DNA forms a prophage). In both cases, the super-infecting or cross-infecting phage DNA is lost, thus contributing to an effective death rate for the phage [we do not include an additional death rate for phage, assuming that they are stable over the timescales of the dynamics of our model (De Paepe and Taddei, 2006)]. Lysogens are known to be immune to super-infections by phages of their type (Echols, 1972) and we expand this into assumption (i) because the two phage variants are identical in all respects except for their propensity to go lytic or lysogenic, i.e., they share immunity proteins. Note that the model doesn't allow formation of double lysogens. Relaxing assumption (ii) doesn't seem to change our conclusions much, even if we allow cross-infections anytime before the decision (see Supplemental Material Section 1.3). In reality, it is quite likely that a cross-infection only affects the decision if it occurs within a short time window after the first infection, in which case the number of cross infected bacteria is likely to be negligible. Finally, we assume that spontaneous induction of lysogens happens very rarely [≤ 10⁻⁵/generation (Bæk et al., 2003)] and would not change the final steady-state much, so we ignore this process.

Maynard et al. (2010) obtained bacterial growth curves for a set of experiments where λ phage infected a population of wild type E. coli and certain mutants. We fitted a version of the above model with just one phage variant (i.e., Equations 1–4 with B₂, L₂ and P₂ set to zero) to several of these growth curves. This is essentially a repeat of what they have done and we find, as they did, that the model fits well (see Supplemental Figures S1, S2). The range of values we observe for each parameter from this set of fits is shown in column 2 of Table 1. We use these values, along with other experimental observations of parameters, to construct a biologically relevant range of values (column 3 of Table 1). Column 4 of the table shows the much larger ranges of parameter values which we explore in our simulations.

We consider the two phage variants in our model as two players in a game, where the winner is defined to be the variant present in larger numbers as lysogens when the dynamics of Equations (1)–(4) reaches steady-state. This is because, in the steady-state, all phage and bacterial populations are zero except for the lysogens, L₁ and L₂, and sometimes the uninfected bacterial population, B₀. We quantify the magnitude of the win by assigning a score or “payoff” of (L₁ − L₂)/(L₁ + L₂) to phage/player 1 (correspondingly, the payoff for phage/player 2 would be (L₂ − L₁)/(L₁ + L₂); notice that the sum of payoffs is always zero, making this, in the language of game theory, a two-player zero-sum game.) To complete the game-theoretic formulation, we need to specify the “strategies” each player can choose. In the simplest case, the strategies correspond to different possible values of the lysogeny propensity. Thus, we use the terms “phage using strategy f” to indicate a phage variant that has lysogeny propensity f. Our model thus consists of a game where we compete two players, phage variant 1 using strategy f₁ and phage variant 2 using strategy f₂, where f₁ and f₂ can take all possible values between 0 and 1. In practice, we calculate the payoffs only for f₁ and f₂ taking the values 0, 0.01, 0.02, …, 0.99, 1, which can be arranged in a “payoff matrix” (see Figure 1A). We later include more complex strategies where the lysogeny propensity is not just a fixed number, but a function of the multiplicity of infection (see next section). Note that, in this paper, we do not explore strategies which would allow, for example, changes of burst size. This is because we want to focus on understanding what determines the lysogeny propensities, so we examine competition between phages that are identical in all respects except for their lysogeny propensity.

Given the payoff matrix, it is possible to define the “best” or “optimal” strategy in different ways. We choose to use the notion of a “minimax” strategy. Figure 1B illustrates how to calculate this from a payoff matrix, but the idea is intuitive: Assuming each player plays the best they can, for each strategy player 2 can play there is a maximum payoff that player 1 can obtain by choosing an appropriate strategy. Therefore, player 2 should play the strategy that minimizes this “maximum payoff” for player 1. In the symmetric zero-sum games we are studying, when there is a unique minimax strategy it has the following property: if both players are playing the minimax strategy and one of them changes its strategy, then the other player can always find a strategy that will beat the player who deviated from the minimax strategy. That is, both players playing the minimax strategy is a sort of “equilibrium”—neither player has an incentive to move away from this.

We extend the model to allow for multiple infections simply by adding more variables that keep track of the population densities of bacteria infected by more than one phage particle (as before, we don't allow cross-infections). We allow multiple infections at any time before the lysis-lysogeny decision. We further assume that decisions happen with the same rate δ irrespective of the multiplicity of infection (MOI).

Here B_i,m is the density of cells infected by m individuals of the ith phage variant. f_i(m) indicates the lysogeny propensity for cells infected by m individuals of the ith variant. In practice, in our simulations, we truncate the multiplicity of infection at 3, which is equivalent to assuming that the propensity of lysogeny is the same for all m ≥ 3.

In this case, the phage strategy is now specified by a set of three lysogeny propensities: f(1) for MOI = 1, f(2) for MOI = 2, and f(3) for MOI ≥3. Here, the number of different strategies is much greater making an exhaustive enumeration of the payoff matrix computationally intensive. Therefore, we instead implement an iterated “evolutionary” game to find the optimal strategy: We begin with two phage variants whose strategies are chosen randomly. That is, for each phage variant, f(1), f(2), and f(3) are independently and randomly chosen from the interval [0 1]. With these values, we run our model consisting of Equations (5)–(9) until the system reaches steady-state, and determine the winner by counting lysogens. In the next iteration of the game, the winner retains its strategy while the loser is replaced by a new player with a “mutated” version of the winner's strategy, where a random change of upto 6% has been made to the lysogeny propensities: f(m) → (1 + 0.01r_m)f(m), where r₁, r₂, and r₃ are random integers uniformly chosen from the interval [−3 3]. Then we run Equations (5)–(9) again with these two phage variants, the winner of the previous iteration vs. its mutant. This procedure of replacing the loser is repeated many times, until the winning strategy stops changing.

The ordinary differential equation (ODE) systems of our phage-bacteria population in Sections 2.4 and 2.5 were simulated in Python v2.7.12 using the native ‘odeint’ solver, which uses Adams or BDF methods to solve non-stiff and stiff ODE systems respectively. We used a custom python script to systematically test all combinations of fixed lysogenic propensity (Section 2.4). The payoff matrix and the minimax value and associated lysogeny propensity were then calculated using MATLAB. The spatial simulations mentioned in Sections 3.2 and Supplemental Material Section 1.6 were implemented with custom code written in C++. The above code files have been deposited in a public Github repository (https://github.com/vaibhhav/phage_competition_paper_frontiers).

We begin by exploring the simplest case where each phage variant is characterized by a fixed propensity of lysogeny. We start the system described by Equations (1)–(4) with a small susceptible bacterial population, at a density B₀(0) well below its carrying capacity, and we introduce the two phage variants in equal, but tiny, amounts P_i(0) [Default starting conditions: P_i(0)=10-4B0(0)]. As mentioned earlier, the idea is to keep all parameters as identical as possible between the two phage variants so that we can focus on the effect of differences in lysogeny propensity alone. For these initial conditions, the steady state typically ends up having only lysogens of both types, with the total bacterial population at carrying capacity (sometimes, e.g., for particularly small burst sizes, uninfected bacteria may also survive in this steady-state; in extreme cases, where one of the lysogeny propensities is zero, there may only be one type of lysogen surviving at the end.)

The dynamics typically goes through three distinct stages (Figure 2):

The buildup: initially, uninfected bacteria grow exponentially, phages grow even faster (but their numbers still remain low enough that they don't affect the growth of uninfected bacteria much), and lysogen numbers are another order of magnitude smaller.

In this steady-state, which phage variant dominates the population can be quantified by the “payoff” (L₁ − L₂)/(L₁ + L₂) for phage variant 1, which ranges from +1 when phage variant 1 completely dominates, to -1 when variant 2 completely dominates (see Section 2.4). For a given parameter set, we run Equations (1)–(4) for all possible combinations of lysogeny propensities f₁ and f₂, and record the payoffs for each case in a matrix. Figure 3 shows this payoff matrix for the parameter values used in Figure 2.

To determine from this the optimal lysogeny propensity, we use a game theoretic view of the phage competition (see Section 2.4), treating each run of Equations (1)–(4) as one two-player game, where the players are the two phage variants and the set of strategies for each player are the different possible values of the lysogeny propensity. From the payoff matrix, we can then calculate the optimal “minimax” strategy where each player aims to play what is best for them, assuming best play by the opponent (see Figure 1). In Figure 3, this optimal strategy lies at the “saddle point” of the payoff matrix (see inset) and corresponds to a lysogeny propensity f_opt ≈ 0.1. For the entire range of parameters we have examined, we find a unique minimax strategy. This means that if both phage variants have lysogeny propensity f_opt, then if either one mutates to have a different lysogeny propensity (whether higher or lower) the other phage variant can always “find” an appropriate lysogeny propensity that will outcompete the phage that mutated away from f_opt—thus neither phage type has any incentive to “choose” a strategy other than f_opt.

The intuitive reason for the existence of a non-zero optimal lysogeny propensity is that aiming to maximize the number of lysogens in the steady-state puts two opposing “forces” on the lysogeny propensity of a phage variant: first, a higher lysogeny propensity increases the number of lysogens formed per infection during the crash phase; second, a lower lysogeny propensity increases the number of phages available for infecting uninfected bacteria during the crash. It is not surprising, therefore, that there is a non-zero lysogeny propensity that balances these opposing forces in our simple ecosystem. Note that the competition is essential for these two opposing forces to exist.

The optimal lysogeny propensity is surprisingly robust to changes in parameter values. Even when the parameters are varied by several orders of magnitude, as shown in Figure 4, we see that f_opt typically lies in the range 5–15%, rising to 30% only when the initial phage population is particularly large. Some trends are visible: f_opt increases with increases in the infection rate constant, η, the burst size, β, the initial phage density, P(0), and the initial bacterial density, B₀(0); and decreases with increases in the decision rate, δ, and the bacterial growth rate, γ (see also Supplemental Figure S4). Allowing cross-infections does not affect this result (see Supplemental Material Section 1.3 and Supplemental Figure S3). Incidentally, even when we relax the well-mixed assumption and run a similar game on a spatial 2d-lattice, we observe similar optimal values of lysogeny propensity; see Supplemental Material Section 1.6 and Supplemental Figure S6.

Using our observation of the three phases of the dynamics, we can explain this remarkable robustness, and the trends, as follows:

In the initial phase, because uninfected bacteria grow exponentially B₀ ~ exp(γt), the free phage abundances grow as a double-exponential with rate dependent on their lytic propensity: log[P_i] ~ (1 − f_i) exp(γt). This assumes that decisions happen quickly, i.e., δ is large compared to γ.

The ratio of lysogens L1/L2=f1P1(t*)/f2P2(t*) will not change from this point until the steady-state is reached, therefore, the lysogeny propensity that maximizes the payoff for phage/player 1 is the value of f₁ that maximizes f1P1(t*) (see Supplemental Material Section 1.4 for a more rigorous analysis that calculates the minimax strategy and gives, in Equation (36) in the supplemental material, essentially the same answer). Using the formula in (a) above, this means:

In other words, the optimal lysogeny propensity depends inversely on t^*, the time of the crash, or equivalently, the duration of the buildup phase. This already explains some of the trends we observe: for instance, a larger burst size or infectivity would cause the crash to happen earlier, and therefore produces a higher f_opt as we observe in Figure 4.

Under the same assumption of a quick decision, we can combine the arguments (a) and (b) above, to calculate how t^* depends on the parameters, finally obtaining (see Supplemental Material Section 1.4):

This explains why the optimal lysogeny propensity is so robust to changes in parameters. Notice that many of the parameters do not even appear in the formula, and the ones that do appear have only a mild effect on the optimal lysogeny propensity because of the logarithm. The burst size is outside the logarithm, but has little effect once it is larger than 10, as it typically is in phages. Supplemental Figure S4 shows that when δ is large, this formula is a good approximation to results obtained from simulations for a large range of parameter values, breaking down only when η, β, B₀(0) and P(0) become very large, or when β and γ are small. When δ is brought down to biologically reasonable values, the formula does not predict f_opt as accurately (the simulations differ by Δf_opt ≈ 0.03) because it ignores the non-negligible time taken to make a decision, but nevertheless it still predicts the robustness and trends well (see Supplemental Figure S4).

Do phages that have a different lysogeny propensity for different multiplicity of infection (MOI) outcompete phages that have a fixed lysogeny propensity that does not depend on MOI? To investigate this, we extend our model to allow multiple infections as described in Section 2.5. A phage “strategy” is now specified by a set of three lysogeny propensities: f(1) for MOI = 1, f(2) for MOI = 2, and f(3) for MOI ≥3 (our model, for simplicity, assumes that the lysogeny propensity is the same for MOI 3 or more). The objective of the games is still the same, to maximize each player's share of lysogens in the steady state. It is computationally rather intensive to exhaustively construct the entire payoff matrix in this case because of the much larger number of possible strategies, so we use a different approach to find the optimal strategy here: Initially, each of the two phage variants starts with random strategies, i.e., randomly chosen values of f(1), f(2), and f(3). With these values, we run our model consisting of Equations (5)–(9) until the system reaches steady-state, and determine the winner by counting lysogens. In the next iteration of the game, the winner retains its strategy while the loser is replaced by a new player with a “mutated” version of the winner's strategy (for details see Section 2.5). These two phage variants now compete again, i.e., we run Equations (5)–(9) once more with these two phage variants, and this procedure of replacing the loser is repeated many times. We find that after many such iterations, the winning strategy, for almost all parameter sets, converges to one where phages always go lytic at MOI 1, and always go lysogenic at MOI 2 or larger. That is, f(1) ≈ 0,while f(2) ≈ f(3) ≈ 1 (see Figure 5). Remarkably, this is very close to what can be inferred from Kourilsky's data: f(1) = 0.005 ± 0.05, f(2) = 0.7 ± 0.3 and f(3) = 0.9 ± 0.1 (Kourilsky, 1973; Avlund et al., 2009b). This was also the optimal strategy we found in exhaustive searches over a restricted strategy space, and we specifically checked that this switch-like strategy outcompetes all phages with fixed lysogeny propensity (see Supplemental Material Section 1.5). Therefore, we believe this strategy is the globally optimal one for most parameter values and initial conditions (very rarely we have observed another switch like strategy, f(1) ≈ f(2) ≈ 0; f(3) ≈ 1, beat this one as shown in Supplemental Figure S5). In the non-well-mixed spatial model we found a similar switch-like strategy was best, except that the switch occurred between MOI = 2 and MOI = 3 (see Supplemental Material Section 1.6 and Supplemental Figure S7).

Our hypothesis for why this MOI-dependent strategy is best amongst the strategies we have examined is that it is able to “detect” the impending crash phase better than other strategies, and respond to it by switching from lysis to lysogeny. Thus, we expect that all other strategies would make more “errors”, i.e., either they would go lysogenic too often in the pre-crash phase or they would go lytic too often during the crash phase. Supplemental Figure S5 provides evidence that this is indeed what typically happens.

Apart from converging to a final state that closely resembles Kourilsky's data, this “evolutionary game,” where the loser is repeatedly replaced by a mutated version of the winner, also highlights certain trends. We observe that, irrespective of the propensity values of the initial random phage strategies, we always see the winning phage's lysogeny propensity for MOI = 1 (i.e., f(1)) rapidly decrease in the first few iterations, while its lysogeny propensity for MOI = 2 (i.e., f(2)) also decreases a little. After f(1) has reached close to zero, in the next iterations, the phage's lysogeny propensity for MOI ≥ 3 (i.e., f(3)) rapidly increases while f(2) also increases, but more slowly. After f(3) has stabilized near 1, f(2) continues to increases till it also reaches 1. In other words, the strongest selection pressure appears to act to make phages always go lytic with single infections. The next strongest selection pressure seems to be for the phages to always go lysogenic when the multiplicity of infection is large. Finally, there is a relatively weaker selection pressure pushing double infections to go predominantly lysogenic. This matches the uncertainty of these lysogeny propensities inferred from Kourilsky's data: his data implies that f(1) is almost certainly very close to 0, and f(3) is very close to 1, while f(2) could really lie anywhere between 0.4 and 1 and still yield a reasonable fit to observations (Kourilsky, 1973; Avlund et al., 2009b).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.