# Evolution of virulence in emerging epidemics: from theory to experimental evolution and back

Wakinyan Benhamou, François Blanquart, Marc Choisy, Thomas Berngruber, Rémi Choquet, Sylvain Gandon

► To cite this version:

Wakinyan Benhamou, François Blanquart, Marc Choisy, Thomas Berngruber, Rémi Choquet, et al. Evolution of virulence in emerging epidemics: from theory to experimental evolution and back. Virus Evolution, 2024, 10 (1), pp.veae069. 10.1093/ve/veae069. hal-04798307

HAL Id: hal-04798307

https://hal.science/hal-04798307v1

Submitted on 9 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Evolution of virulence in emerging epidemics: from theory to experimental evolution and back

Wakinyan Benhamou<sup>1,†,*</sup>, François Blanquart<sup>2,†</sup>, Marc Choisy<sup>3,4</sup>, Thomas W. Berngruber<sup>5</sup>, Rémi Choquet<sup>1,‡</sup>, Sylvain Gandon<sup>1,‡</sup>

## Abstract

The experimental validation of theoretical predictions is a crucial step in demonstrating the predictive power of a model. While quantitative validations are common in infectious diseases epidemiology, experimental microbiology primarily focuses on the evaluation of a qualitative match between model predictions and experiments. In this study, we develop a method to deepen the quantitative validation process with a polymorphic viral population. We analyse the data from an experiment carried out to monitor the evolution of the temperate bacteriophage λ spreading in continuous cultures of *Escherichia coli*. This experimental work confirmed the influence of the epidemiological dynamics on the evolution of transmission and virulence of the virus. A variant with larger propensity to lyse bacterial cells was favoured in emerging epidemics (when the density of susceptible cells was large), but counter-selected when most cells were infected. Although this approach qualitatively validated an important theoretical prediction, no attempt was made to fit the model to the data nor to further develop the model to improve the goodness of fit. Here, we show how theoretical analysis—including calculations of the selection gradients—and model fitting can be used to estimate key parameters of the phage life cycle and yield new insights on the evolutionary epidemiology of the phage λ. First, we show that modelling explicitly the infected bacterial cells which will eventually be lysed improves the fit of the transient dynamics of the model to the data. Second, we carry out a theoretical analysis that yields useful approximations that capture at the onset and at the end of an epidemic the effects of epidemiological dynamics on selection and differentiation across distinct life stages of the virus. Finally, we estimate key phenotypic traits characterizing the two strains of the virus used in our experiment such as the rates of prophage reactivation or the probabilities of lysogenization. This study illustrates the synergy between experimental, theoretical, and statistical approaches; and especially how interpreting the temporal variation in the selection gradient and the differentiation across distinct life stages of a novel variant is a powerful tool to elucidate the evolutionary epidemiology of emerging infectious diseases.

**Keywords:** evolutionary epidemiology; life-history evolution; competition; natural selection; statistical inference; bacteriophage λ

## 1. Introduction

Evolutionary epidemiology theory predicts that the evolution of pathogen transmission is driven by the availability of susceptible hosts. At the onset of an epidemic, when the density of susceptible hosts is high, more transmission is favoured by natural selection. When a positive covariance exists between transmission and virulence (Anderson and May, 1982; Alizon et al., 2009; Alizon and Michalakis, 2015), this selection for higher transmission can indirectly select for higher virulence (Bull, 1994; Day, 2002; Lenski and May, 1994; Frank, 1996; Day and Proulx, 2004; Gandon and Day, 2007). Yet, an experimental validation of this prediction was needed to demonstrate the relevance of these predictions on the evolution of pathogens in emerging epidemics.

This prediction was put to the test in a previous study using experimental evolution of the temperate bacteriophage (or phage) λ (Berngruber et al., 2013). Phages are viruses that infect bacteria and phage λ is the archetypal temperate phage, which can switch between a lytic and a lysogenic life style. Upon infection, the virus may commit to the lytic pathway by hijacking the host's replication machinery to produce new virions (viral particles) and eventually release them in the environment after the lysis of the host cell. Alternatively, the virus may commit to the lysogenic pathway by integrating its genome into the bacterial chromosome where it will lie in a dormant state as a prophage and be hereditarily transmitted to daughter cells at the pace of lysogen divisions. The prophage may also regain virulence by excising itself from the

<footnote n="1"><sup>1</sup>CEFE, Univ Montpellier, CNRS, EPHE, IRD, Montpellier, France</footnote>

<footnote n="2"><sup>2</sup>Centre for Interdisciplinary Research in Biology (CIRB), Collège de France, CNRS, INSERM, PSL Research University, Paris, France</footnote>

<footnote n="3"><sup>3</sup>Centre for Tropical Medicine and Global Health, Nuffield Department of Medicine, University of Oxford, Oxford OX3 7BN, United Kingdom</footnote>

<footnote n="4"><sup>4</sup>Oxford University Clinical Research Unit, Ho Chi Minh City, Vietnam</footnote>

<footnote n="5"><sup>5</sup>Reinhard-zu-Rhynernstrasse 7, Hamm 59069, Germany</footnote>

<footnote><sup>†</sup>These authors contributed equally to this work.</footnote>

<footnote><sup>‡</sup>These authors also contributed equally to this work.</footnote>

<footnote>*Corresponding author. CEFE, Univ Montpellier, CNRS, EPHE, IRD, Montpellier, France. E-mail: wakinyan.benhamou@cefe.cnrs.fr</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Flow chart of the phage-bacteria system. The subscripts *w* and *m* refer to the wildtype (in green) and mutant (or virulent, in red) strain of phage λ, respectively. In the bacterial population (*E. coli*) of size *N*, bacteria are either susceptible (*S*), lysogenic (*L*)—i.e. carrying a prophage in their chromosome (small ellipses)—or carrying phages (small pentagons) replicating in their cytoplasm prior to lysis (*Y*); *V* is the free virus stage (culture medium). Width of the arrows between *S*, *L*, and *Y* reflects the relative rates of different events: the wildtype strain is transmitted mostly vertically (high probability of lysogenization and low rate of prophage reactivation) while the virulent strain is transmitted mostly horizontally (low probability of lysogenization and high rate of prophage reactivation). Dashed arrows symbolize the role of the free viral particles in the force of infection (epidemiological feedback). See notations in Table 1.</figcaption></figure>

host genome, switching to a lytic cycle (reactivation, also called
induction) and thus shifting from vertical to horizontal transmis-
sion (Lwoff, 1953; Echols, 1972; Ptashne, 1992; Gandon, 2016). The
evolution experiment designed in Berngruber et al. (2013) moni-
tored the competition between two strains of phage λ with distinct
life-history strategies in continuous cultures of Escherichia coli. The
first strain is the wildtype, which is known to have a relatively large
lysogenization rate and low reactivation rate. The second strain is
the λcl857 variant, which carries a point mutation in the transcrip-
tional repressor protein cl (St-Pierre and Endy, 2008; Sussman and
Jacob, 1962); the cl mutant is known to be more virulent and trans-
mitted mostly horizontally through lytic cycles. Berngruber et al.
(2013) developed a mathematical model tailored to the life cycle
of phage λ. Numerical simulations of this model using parameter
values from previous experimental studies led to three theoretical
predictions: (i) the virulent strain outcompetes the wildtype when
susceptible hosts are abundant, but the direction of selection is
reversed as soon as the epidemic reaches high prevalence, (ii) the
lower the initial prevalence, the higher the increase in virulence
during the epidemic and (iii) the virulent strain is always more
frequent among viral particles than among prophages. Tracking
both the epidemiology (prevalence) and the evolution of the virus
(frequency of the virulent strain among viral particles and among
infected bacteria), all three predictions were confirmed experi-
mentally (Berngruber et al., 2013). Yet, the data were only used
as a qualitative validation of the theory and no attempt was made
to explore the quantitative match between the predicted and the
observed dynamics of the virus.

In the present work, we show how the quantitative analysis of
the experimental results from Berngruber et al. (2013) improves
our understanding of the evolutionary epidemiology of phage λ.
First, we modified the structure of the epidemiological model to
better capture the transient evolutionary dynamics of the virus
among infected bacteria. Second, we carry out a theoretical anal-
ysis of this model to provide useful approximations to predict
the evolutionary dynamics of the virus at different stages of the
epidemic. In particular, we compute the selection gradients and

the differentiation across distinct life stages of the virus. Finally,
we develop a statistical inference approach to obtain quantita-
tive estimates of the parameters of the model and, especially, the
life-history traits of the different strains of phage λ.

## 2. Materials and methods

Following Berngruber et al. (2013), we first model the competition
between the wildtype strain—hereafter denoted by the subscript
*w*—and the mutant strain λcl857—hereafter denoted by the sub-
script *m*—of phage λ in a chemostat with a well-mixed continuous
culture of its bacterial host *E. coli* (Fig. 1). We summarize the nota-
tions in Table 1. We then recall the experimental data generated in
the original study and describe how we generate simulated data
to validate our ability to estimate model parameters. In the last
section, we detail the statistical inference approach we developed.

## 2.1 A model coupling epidemiology and evolution of phage virulence

### 2.1.1 Epidemiology

Bacteria are either susceptible (S) or infected with phage λ (we
do not consider resistant bacteria); infected bacteria may be lyso-
genic (L), following phage integration (lysogenization), or carrying
phages replicating in their cytoplasm prior to lysis (Y), following
either lytic infection or prophage reactivation. Unlike the original
model in Berngruber et al. (2013), adding an extra stage Y allows us
to take the lysis time into account; we show below that this new
model improves the goodness of fit to the data (see Results section
'3.2.2 Inference from experimental data'). The free virus stage (V)
corresponds to the free viral particles (virions) that are in the cul-
ture medium—'free' meaning 'extracellular' here. Throughout, for
each state variable (aka compartment), for instance S, we denote
by S(t) its density at the current time t and Š(t) its derivative with
respect to time. Hence, L(t) = Lw(t) + Lm(t) is the total density at time
t of L cells, Y(t) = Yw(t) + Ym(t), of Y cells, V(t) = Vw(t) + Vm(t), of free
viral particles, and N(t) = S(t) + L(t) + Y(t), of bacteria.
<table><caption>**Table 1.** Notations. The subscripts *w* and *m* refer to the wild-type strain and the mutant (or virulent) strain λcl857 of phage λ, respectively. Overlines refer to mean values of life-history traits across all genotypes. SD stands for 'standard deviation'.</caption><thead><tr><th>Term</th><th>Definition</th></tr></thead><tbody><tr><td>N</td><td>Bacterial population</td></tr><tr><td>S</td><td>Susceptible bacteria</td></tr><tr><td>L, L<sub>w</sub>, L<sub>m</sub></td><td>Lysogens (bacteria carrying a prophage)</td></tr><tr><td>Y, Y<sub>w</sub>, Y<sub>m</sub></td><td>Bacteria with phages replicating in their cytoplasm prior to lysis</td></tr><tr><td>V, V<sub>w</sub>, V<sub>m</sub></td><td>Free viral particles (virions)</td></tr><tr><td>P</td><td>Prevalence (1-S/N)</td></tr><tr><td>p</td><td>Frequency of L cells infected by the mutant strain m</td></tr><tr><td>q</td><td>Frequency of the mutant strain m in the free virus stage (V)</td></tr><tr><td>f</td><td>Frequency of Y cells infected by the mutant strain m</td></tr><tr><td>g</td><td>Frequency of cells infected (either Y or L) by the mutant strain m</td></tr><tr><td>α&#773;, α<sub>w</sub>, α<sub>m</sub></td><td>Rates of prophage reactivation; α<sub>w</sub> &lt; α<sub>m</sub></td></tr><tr><td>φ&#773;, φ<sub>w</sub>, φ<sub>m</sub></td><td>Probabilities of lysogenization (phage integration into the host genome upon infection); φ<sub>w</sub> &gt; φ<sub>m</sub></td></tr><tr><td>Δα, Δφ</td><td>Phenotypic differences between the mutant and wildtype strain; Δα = α<sub>m</sub> - α<sub>w</sub>, Δφ = φ<sub>m</sub> - φ<sub>w</sub></td></tr><tr><td>a</td><td>Adsorption rate of free viral particles onto the surface of bacteria</td></tr><tr><td>b</td><td>Probability of fusion (injection of the phage genome) upon adsorption</td></tr><tr><td>τ</td><td>Lysis rate; 1/τ corresponds to the mean lysis time</td></tr><tr><td>B</td><td>Burst size (number of released viral particles upon lysis)</td></tr><tr><td>r</td><td>Bacterial intrinsic growth rate</td></tr><tr><td>K</td><td>Carrying capacity for the bacterial population: K = 10<sup>9</sup> cells</td></tr><tr><td>δ</td><td>Dilution rate of the continuous culture: δ = 0.8 h<sup>-1</sup></td></tr><tr><td>R<sub>0,w</sub>, R<sub>0,m</sub></td><td>Basic reproduction numbers of the phage</td></tr><tr><td>S</td><td>Selection gradient of the virulent phage (rate at which it grows or declines in frequency on the logit scale)</td></tr><tr><td>Q<sup>VL</sup></td><td>Differentiation of the virulent phage between free phages and prophages</td></tr><tr><td>σ<sub>P</sub>, σ<sub>g</sub>, σ<sub>q</sub></td><td>SD of measurement errors for logit(P), logit(g) and logit(q), respectively</td></tr></tbody></table>

Susceptible and lysogenic bacteria grow at a *per capita* logistic rate $r(1-N(t)/K)$, where *r* is the intrinsic growth rate and *K* the carrying capacity. We assume, as in *Berngruber et al.* (2013), that the prophage does not affect the intrinsic growth rate of its host and that vertical transmission is perfect. Bacteria and virions are removed from the chemostat at a dilution rate $\delta$. Free viral particles adsorb onto bacterial cells at a rate $\alpha$; adsorption is non-reversible, that is, the fate of adsorbed viruses is only to infect or die. Infection of a susceptible host also requires the injection of the phage's genetic material into the bacterial cytoplasm (fusion, with probability *b*). The force of infection—i.e. the *per capita* infection rate—is therefore given by $abV(t)$. We assume that super-infection (including coinfection with both strains) does not occur. In particular, prophage establishment of phage $\lambda$ is known to provide cellular immunity, or superinfection inhibition (Lwoff, 1953; Ptashne, 1992; *Berngruber et al.*, 2010; Gandon, 2016). Upon infection, the wildtype and virulent strains of phage $\lambda$ may either be integrated as prophages into the host genome (lysogenic cycle) with probabilities $\phi_w$ and $\phi_m$, respectively, such that $\Delta\phi = \phi_m - \phi_w < 0$, or start the biosynthesis and assembly of viral copies (lytic cycle) with complementary probabilities. Once integrated, reactivations of the wildtype and virulent prophages occur at a rate $\alpha_w$ and $\alpha_m$, respectively, such that $\Delta\alpha = \alpha_m - \alpha_w > 0$. Following either lytic infections or prophage reactivations, host cells are lysed at a rate $\tau$ and eventually release B viral particles upon lysis (burst size). We

assume that parameters *a*, *b*, $\tau$, and *B* are the same for the wild-type and the mutant which yields the following system of ordinary differential equations (ODEs):

$$ 
\left\{
\begin{aligned}
\dot{S}(t) &= rS(t)\left(1 - \frac{N(t)}{K}\right) - \underbrace{abV(t)S(t)}_{\text{Growth}} - \underbrace{\delta S(t)}_{\text{Removal}} \\
\dot{L}(t) &= rL(t)\left(1 - \frac{N(t)}{K}\right) + \underbrace{\phi(t)abV(t)S(t)}_{\text{Growth & vertical transmission}} - \underbrace{\bar{\alpha}(t)L(t)}_{\text{Prophage reactivation}} - \underbrace{\delta L(t)}_{\text{Removal}} \\
\dot{Y}(t) &= \underbrace{(1 - \phi(t))abV(t)S(t)}_{\text{Lytic infection}} + \underbrace{\bar{\alpha}(t)L(t)}_{\text{Prophage reactivation}} - \tau\dot{Y}(t) - \underbrace{\delta Y(t)}_{\text{Lysis}} \\
\dot{V}(t) &= \tau\dot{Y}(t)\mathbf{B} - \underbrace{aN(t)V(t)}_{\text{Virion release}} - \underbrace{\delta V(t)}_{\text{Adsorption Removal}}
\end{aligned}
\right.
, \quad (1) $$

where $\bar{\phi}(t)$ is the mean probability of lysogenization upon infection and $\bar{\alpha}(t)$ the mean rate of reactivation among prophages:

$$ \left\{ \begin{aligned} \bar{\phi}(t) &= q(t)\phi_m + (1-q(t))\phi_w, \\ \bar{\alpha}(t) &= p(t)\alpha_m + (1-p(t))\alpha_w, \end{aligned} \right. \quad (2) $$

with $q(t) = V_m(t)/V(t)$ and $p(t) = L_m(t)/L(t)$, the frequencies of the virulent strain at time *t* among compartments V and L, respectively.

### 2.1.2 Evolution

Along with these epidemiological dynamics, we also track the evolutionary dynamics of phage $\lambda$. We recall that $p(t)$ and $q(t)$ refer to frequencies of the virulent strain in compartment L and V, respectively. In addition, we also denote $f(t) = Y_m(t)/Y(t)$, the frequency of the virulent strain at time *t* among Y cells, and $g(t)$, the frequency of the virulent strain in infected cells (either Y or L) such that:

$$ g(t) = \frac{Y_m(t) + L_m(t)}{Y(t) + L(t)} = f(t) \left( \frac{Y(t)}{Y(t) + L(t)} \right) + p(t) \left( \frac{L(t)}{Y(t) + L(t)} \right). \quad (3) $$

Using (1) and (2), one may calculate the ODE for each of these frequencies (see Supplementary Appendix S1.2). More conveniently, we will then focus on logit-frequencies instead, that is, the log odds ln(frequency of the mutant/frequency of the wild-type). Taken together, the equations of the temporal dynamics of the frequencies and the model (1) yield the coupled evolutionary-epidemiological dynamics of this phage-bacteria system. The analysis of this model can provide key insight on the evolutionary forces acting on the virus; it may also provide useful approximations for the change in mutant frequency at different stages of the epidemic.

## 2.2 Time series datasets

### 2.2.1 Experimental data

We use experimental time series obtained in the first evolution experiment of *Berngruber et al.* (2013). Briefly, this experiment started with eight independent chemostats (5 mL chamber volume, maintained at 35°C) of well-mixed bacterial cultures of E. coli MG1655 (RecA+) at carrying capacity. Bacteria were initially infected by both the wildtype and virulent strain λcl857 of phage λ (prophage stage with initial ratio 1:1). Two treatments were considered using four chemostats each: (i) an epidemic treatment—low initial prevalence, around 1%—and (ii) an endemic treatment—high initial prevalence, around 99%. Throughout the course of the experiment, several quantities were monitored: the prevalence, the frequency infected by each strain and the strain frequency in the culture medium (free virus stage). Samplings in each chemostat were performed hourly, from 1 to 60 h maximum. The prevalence and the frequency of infected hosts were tracked using flow
cytometry (FACS) with fluorescent protein marker colours (CFP and YFP) while the strain frequency in the free virus stage was tracked by qPCR.

## 2.2.2 Simulated data

Alongside experimental data, we also carry out an analysis based on simulated data in order to validate our ability to infer parameters from experimental data. For this purpose, we take: $\alpha_w = 7 \times 10^{-3}$, $\alpha_m = 2 \times 10^{-2}$, $\phi_w = 0.2$, $\phi_m = 2 \times 10^{-2}$, $a = 3 \times 10^{-9}$, $b = 0.1$, $B=80$, $r = 1.4$, $\tau = 1.5$, $K = 10^9$ and $\delta = 0.8$. At $t=0$, bacteria are at carrying capacity K with initial prevalence 1% (epidemic treatment) or 99% (endemic treatment) and the initial prophage ratio for the two strains is 1:1. We simulate the deterministic model (1) for the epidemic and endemic treatment (Supplementary Fig. S1). We then add i.i.d. Gaussian noise at each time point to mimic measurement errors on the logit-prevalence logit($P(t)$), the logit-frequency of hosts infected by the virulent phage logit($g(t)$) and the logit-frequency of the virulent phage in the free virus stage logit($q(t)$)—this is independently repeated four times for each simulation to obtain four replicates (chemostats) per treatment. We modulate data quantity through two sampling frequencies: $0.1 h^{-1}$ vs. $1 h^{-1}$, and we modulate data quality through two standard deviations (SD) of measurement errors: $0.01$ vs. $0.5$. We thus end up with four combinations of data quantity and quality (see example in Supplementary Fig. S2) from which we then try to recover parameter values.

# 2.3 Maximum likelihood estimation

For the estimation process, we used a two-step approach. Using theoretical results, we first compute point estimates of the rates of prophage reactivation $\alpha_w$ and $\alpha_m$. Then, we fix the latter to estimate the remaining parameters of the model using non-linear optimizations. Note that it is possible to run non-linear optimizations to estimate all parameters but this two-step approach makes the optimization easier by reducing the dimensionality of the problem. The Sieve bootstrap method (Bühlmann, 1997; Ulloa et al., 2013) is used to compute the joint distributions of all these estimated parameters.

## 2.3.1 Estimation of the rates of prophage reactivation

From the analysis of the model (see Results section '3.1 Theoretical analysis'), we show that when the system reaches high prevalence the selection gradient $S$ of the mutant—i.e. the rate at which it grows or declines in frequency on the logit scale—is simply given by: $S = \alpha_w - \alpha_m = -\Delta\alpha$ (see Results section '3.1.2 Evolution at the end of the epidemic'). Furthermore, the differentiation $Q^{VL}$ of the virulent strain between free phages and prophages converges towards approximately: $Q^{VL} = \alpha_m/\alpha_w = 1 + \Delta\alpha/\alpha_w$ (see Results section '3.1.3 Differentiation across compartments'). Combining these two expressions enables us to estimate separately both reactivation rates:

$$ \begin{cases} \alpha_w = \frac{\delta}{1 - Q^{VL}} \\ \alpha_m = \frac{\delta \times Q^{VL}}{1 - Q^{VL}}. \end{cases} $$

For each chemostat, we therefore only keep the data from the time point the prevalence has reached 95%. We fit a linear model on the logit-frequency infected by the virulent phage logit($g(t)$) to estimate the slope $S$, and we estimate $Q^{VL}$ by calculating the geometric mean of $\frac{q(t)(1-q(t))}{(1-q(t))g(t)}$ (see details in Supplementary Appendix S3.1). Note that we substitute $p(t)$ by $g(t)$ for the calculation of the

differentiation because we only have access to $g(t)$ in the experiment, and $g(t)$ is almost identical to $p(t)$ towards the end of the epidemic (see Results section '3.1.2 Evolution at the end of the epidemic' and Supplementary Fig. S1).

## 2.3.2 Estimation of the remaining parameters

We have three response variables: (i) logit($P(t)$), the logit-prevalence, (ii) logit($g(t)$), the logit-frequency of hosts infected by the virulent phage, and (iii) logit($q(t)$), the logit-frequency of the virulent phage in the free virus stage. Let $P_0^{epidemic}$ and $P_0^{endemic}$ be the initial conditions (at $t=0$) of the prevalence in the epidemic and endemic treatments, respectively, and $p_0$, the initial condition of the frequency infected by the virulent phage (identical for both treatments). Let $\theta = (\phi_w, \phi_m, a, b, r, \tau, P_0^{epidemic}, P_0^{endemic}, p_0)$ be the vector of model parameters to estimate. Several parameter values are fixed: $K = 10^9$ cells and $\delta = 0.8 h^{-1}$ (Berngruber et al., 2013); $\alpha_w$ and $\alpha_m$ are fixed to their previous point estimates (cf. Material and methods section '2.3.1 Estimation of the rates of prophage reactivation'); and $B = 80$ virus cell$^{-1}$ (Wang, 2006) as we show in Results section '3.2.1 Inference from simulated data' that the burst size B is not separately identifiable from parameter b. We estimate as well $\sigma_p$, $\sigma_g$, and $\sigma_q$, the SD of measurement errors for each of the three response variables, respectively. We assume that variation around the deterministic dynamics stems solely from measurement errors and we thereby neglect any additional process stochasticity. This is justified in particular by the very controlled conditions of the experiment and by the large population sizes of both the bacteria and the phage in the chemostats. For each response variable, measurement errors are assumed to be normally distributed and i.i.d. across all treatments, replicates, and time points. The corresponding likelihoods ($\mathcal{L}$) are thence respectively given by:

$$ \mathcal{L}_P(\theta, \sigma_p) = \prod_{i,j,t} \varphi \left( \logit(P_{ij}(t))^{\text{data}} \mid \logit(P_{ij}(\theta, t))^{\text{sim}}, \sigma_p^2 \right), $$

$$ \mathcal{L}_g(\theta, \sigma_g) = \prod_{i,j,t} \varphi \left( \logit(g_{ij}(t))^{\text{data}} \mid \logit(g_{ij}(\theta, t))^{\text{sim}}, \sigma_g^2 \right), $$

$$ \mathcal{L}_q(\theta, \sigma_q) = \prod_{i,j,t} \varphi \left( \logit(q_{ij}(t))^{\text{data}} \mid \logit(q_{ij}(\theta, t))^{\text{sim}}, \sigma_q^2 \right), $$

where i refers to treatments (epidemic vs. endemic), j to the jth replicate (chemostat), t to time points and $\varphi(\mid \mu, \sigma^2)$ to the probability density function of the Normal distribution with mean $\mu$ and variance $\sigma^2$; sim indicates model outputs. We then denote $\hat{\theta}$, the maximum likelihood estimation (MLE) estimator of $\theta$, such that:

$$ \hat{\theta} = \underset{\theta}{\operatorname{argmax}} \left( \ln(\mathcal{L}_P(\theta, \sigma_p)) + \ln(\mathcal{L}_g(\theta, \sigma_g)) + \ln(\mathcal{L}_q(\theta, \sigma_q)) \right). $$

In practice, we minimize the negative overall log-likelihood using the Nelder-Mead (aka downhill simplex) algorithm (Nelder and Mead, 1965). Parameter bounds (reported in Supplementary Table S2) are enforced through parameter transformations. Due to the presence of local minima, optimizations are repeated for 2000 sets of uniformly drawn starting points to ensure convergence to a global minimum. The best MLE set of estimates $\hat{\theta}$ corresponds to the fit associated with the lowest negative overall log-likelihood with successful completion.

## 2.3.3 Confidence intervals

We generate bootstrapped data to compute 95% CIs of our parameters using Sieve bootstrap (Bühlmann, 1997; Ulloa et al., 2013) on the residuals between experimental data and the best fit
of our model. For this purpose, autoregressive moving-average (ARMA) models are fitted to the time series of centred residuals of each chemostat independently. We then use ARMA models to simulate new residuals from which we reconstruct new datasets. We eventually reiterate the above estimation procedure (cf. Materials and methods section '2.3.1 Estimation of the rates of prophage reactivation' to '2.3.2 Estimation of the remaining parameters'), but starting non-linear optimizations only from the best MLE estimates we obtained with the original data. By repeating this for 10 000 bootstrapped datasets, we compute the joint distributions of estimated parameters.

## 2.4 Details of the implementation

Numerical simulations and data analyses were carried out using R (R Core Team, 2022) version 4.2.0 (2022-04-22). ODEs were solved numerically by the function ode—with method lsoda—from the package deSolve (Soetaert et al., 2010). Non-linear optimizations for MLE were tackled with the function nmk (Kelley, 1999), from the package dfoptim, which gave here more stable results than optim from base R. Fit and selection of ARMA models for Sieve bootstrap were carried out using the function auto.arima, from the package forecast (Hyndman and Khandakar, 2008).

# 3. Results

## 3.1 Theoretical analysis

When $r > \delta$, the virus-free system converges to an equilibrium where $S(\infty) = K(1 - \delta/r)$. When a single strain $k \in \{w, m\}$ of the virus (with phenotypes $\alpha_k$ and $\phi_k$) is introduced in the bacterial population (fully susceptible, with density $S_0$), the fate of the phage-bacteria system depends on the basic reproduction number of the pathogen $\mathcal{R}_{0,k}$ which is given by:

$$ \mathcal{R}_{0,k} = \frac{A + \sqrt{A^2 - 4r \left(1 - \frac{S_0}{K}\right) (1 - \phi_k) abS_0 \tau B (aS_0 + \delta)(\alpha_k + \delta)(\tau + \delta)}}{2(aS_0 + \delta)(\alpha_k + \delta)(\tau + \delta)}, \quad (4) $$

with $A = r(1 - \frac{S_0}{K})(aS_0 + \delta)(\tau + \delta) + abS_0\tau B(\alpha_k + (1 - \phi_k)\delta)$ [see Supplementary Appendix S2.1 for the construction of the next-generation matrix, following Diekmann et al. (2010)]. When $\mathcal{R}_{0,k} < 1$, the virus goes extinct and the bacterial population converges to the previous virus-free equilibrium. Alternatively, when $\mathcal{R}_{0,k} > 1$, an epidemic breaks out and eventually stabilizes to an endemic equilibrium where all the cells are infected by the virus (Supplementary Appendix S2.1).

In the following, we analyse the evolutionary dynamics (i) at the beginning of an epidemic where $S(t) = S_0$ and (ii) when the system stabilizes towards the endemic equilibrium where $S(t) = 0$.

### 3.1.1 Evolution in an emerging epidemic

At the onset of the epidemic, susceptible cells are highly abundant. For the sake of simplicity, we analyse the dynamics of the epidemic when the host density is assumed to be constant over time ($S(t) = S_0$). As in the experimental design, we start with a bacterial population at carrying capacity ($S_0 \approx K$). Density dependence reduces cell reproduction and vertical transmission of the virus. Consequently, the epidemic is mainly driven by the lytic pathway. At $t=0$ though, only lysogens are introduced at very low density. After a very short time, the phage-bacteria system reaches its new dynamical regime, following prophage reactivations, lyses and releases of virions. From there, we assume that $L(t)/Y(t) \approx 0$ and we focus on the dynamics of the mutant in the compartment

Y. We derive an approximation of the selection gradient $\mathcal{S}$ of the virulent phage, which corresponds to the rate at which it increases in frequency on the logit scale, and we show in Supplementary Appendix S2.2 that:

$$ \mathcal{S} \propto -\Delta\phi, \quad (5) $$

meaning that it is selected for ($\Delta\phi < 0$, cf. Supplementary Fig. S3A). This is consistent with the experimental data where the virulent phage transiently outcompetes the wildtype at the early stage of the epidemic treatment (Berngruber et al., 2013). Our prediction is accurate when the density of susceptible hosts remains effectively constant over time; otherwise, however, our prediction deviates from the simulation whose rate slows down because the density of susceptible cells is rapidly decreased by the spread of the epidemic (Fig. 2A and Supplementary Fig. S3C).

### 3.1.2 Evolution at the end of the epidemic

At the end of the epidemic, we expect that all the cells will be infected by a prophage and there will no longer be any susceptible cells. Consequently, no horizontal transmission takes place and, in contrast with the previous scenario, we can neglect the density of Y cells relative to the density of lysogens. Indeed, since $1/\alpha_w \gg 1/\tau$ (time elapsed between phage integration and reactivation is much longer than lysis time), then $Y(t)/L(t) \approx 0$ (see details in Supplementary Appendix S2.3). Note that this also means that the frequency infected by the virulent strain is mainly driven by the corresponding frequency in L, that is, $g(t) \approx p(t)$. We show in Supplementary Appendix S2.3 that the selection gradient of the virulent phage is given by:

$$ \mathcal{S} = -\Delta\alpha. \quad (6) $$

As a result, the virulent phage is counter-selected in the long-term ($\Delta\alpha > 0$) and, in each compartment, linearly declines in frequency at a rate $\mathcal{S}$ (negative slope) on the logit scale (Fig. 2B and Supplementary Fig. S1B).

### 3.1.3 Differentiation across compartments

The two previous subsections deal with the dynamics of evolution over time; let us now look at the dynamics of selection between different compartments. We define here the differentiation $Q^{VL}$ between the free virus stage and lysogens as:

$$ Q^{VL}(t) = \frac{q(t)}{1 - q(t)} \frac{1 - p(t)}{p(t)}, \quad (7) $$

such that:

$$ \ln(Q^{VL}(t)) = \logit(q(t)) - \logit(p(t)). $$

We show in Supplementary Appendix S2.4 that, at the end of the epidemic, the differentiation $Q^{VL}(t)$ converges approximately towards:

$$ Q^{VL} \approx \frac{\alpha_m}{\alpha_w} = 1 + \frac{\Delta\alpha}{\alpha_w}, \quad (8) $$

(Fig. 2C), a value of 1 corresponding to no differentiation. We obtain the same expression for the convergence of the differentiation between Y and L cells and around 1 between the free virus stage V and Y cells (cf. Supplementary Appendix 2.4). Interestingly, these results are also valid at $t = 0^+$ (Supplementary Fig. S4) as there are no free phage particles yet and thereby no horizontal transmission which, within a very short space of time, resembles the end of the epidemic. When the system stabilizes at the end of
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Theoretical predictions. At $t=0$, bacteria are at carrying capacity K with initial prevalence 1% (epidemic treatment) or 99% (endemic treatment). The initial phage ratio for the two strains is 1:1. See Supplementary Table S1 for parameter values. (a) The virulent phage is selected for at the early stage of the epidemic. Our approximation for the trajectory of the logit-frequency of the virulent phage (light red shading, see Supplementary Appendix S2.2) predicts the dynamics of the virulent phage when the density of susceptible cells is kept constant (dashed red line). The slope of the upper bound of our approximation is proportional to $-\Delta\phi = \phi_w - \phi_m$ (equation (5)). As only lysogens are introduced in the susceptible bacterial population at $t=0$, we let the phage-bacteria system reach its new dynamical regime and only start our approximation at $t=1$ (white dot). However, the rapid drop in the availability of susceptible hosts due to the viral epidemic reduces the increase of the virulent phage relative to our approximation (compare the full red line and the dashed red line). (b) The virulent phage is counter-selected when the epidemic reaches high prevalence (indicated with the white dots in the epidemic and endemic treatments). At the end of the epidemic, the logit-frequency of the virulent phages decreases linearly with negative slope $-\Delta\alpha = \alpha_w - \alpha_m$ (equation (6)). (c) The virulent strain is more frequent among free viruses (V) than among prophages (L). When the system reaches high prevalence, we predict the differentiation between these two compartments to converge towards an equilibrium that is approximately equal to $1 + \Delta\alpha/\alpha_w = \alpha_m/\alpha_w$ (equation (8)).</figcaption></figure>

the epidemic, the virulent strain is therefore more frequent among free viruses and Y cells than among L cells (prophages) but we also expect almost no differentiation between free viruses and Y cells (Supplementary Fig. S4).

## 3.2 Statistical inference

### 3.2.1 Inference from simulated data

We first conduct a simulation study. We start by validating our estimation method of the rates of prophage reactivation $\alpha_w$ and

$\alpha_m$. Combining equations (6) and (8), we compute point estimates of both parameters (cf. Materials and methods section '2.3.1 Estimation of the rates of prophage reactivation') for a large number of simulated datasets. Overall, estimated values show a very good match with those used in the original simulations (Supplementary Fig. S5), especially when the SD of measurement errors is small.

We then show that parameters *b* and *B* may not be separately identifiable using the simulated dataset closest to the original deterministic simulation (sampling effort = 10 h⁻¹ and SD of measurement error = 0.01, see Supplementary Fig. S2A). We fix K
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Differentiation of the virulent strain between free phages and infected cells. We compute $q(t)(1-q(t))/(1-q(t)g(t))$ for each chemostat (thin lines) and the corresponding geometric mean for each treatment (thick lines). The value 1 (horizontal grey line) corresponds to no differentiation. For each chemostat, values before the time point the system has reached a prevalence of 95% are shown in transparent color. With high prevalence, we have $g(t) \approx p(t)$ and we thereby obtain the differentiation of the virulent strain between free phage particles (V) and prophages (L). In these conditions, such differentiation is expected to reach an equilibrium—the horizontal dashed line indicates the corresponding geometric mean (4.60) across all chemostats and time points—with theoretical value $1 + \Delta\alpha/\alpha_w = \alpha_m/\alpha_w$ (cf. Fig. 2C).</figcaption></figure>

and $\delta$ to their known values, as well as $\alpha_w$ and $\alpha_m$, as though the rates of prophage reactivation have been correctly estimated beforehand. Setting (b, B) to many combinations of values ((b, B) $\in$ [0, 0.2] $\times$ [0, 100]) and maximizing over the remaining parameters, the resulting landscape of the overall log-likelihood suggests that only the product $b \times B$ is identifiable (Supplementary Fig. S6). Such compensatory effect may be understood as the impossibility here of discriminating whether phage $\lambda$ infect fewer bacteria (lower probability $b$ of fusion) but release more viral copies upon lysis (larger burst size B) or vice-versa. This is the reason why we fix B in the following and only estimate b.

Fixing B to its true value used in the original simulations, point estimates of parameters are quite close to their true values (Supplementary Fig. S7). However, some parameters—like $r$, $\tau$ or $b$—seem more difficult to estimate with a lower sampling effort.

### 3.2.2 Inference from experimental data

The original model used in Berngruber et al. (2013) failed to properly capture the evolutionary dynamics among infected hosts at the early stage of the epidemic. Such discrepancy was an opportunity to go back from experimental data to theory and to challenge the structure of the model. We noticed that adding an extra stage Y in the infected compartment allowed us to better fit the experimental data than the original model ($\Delta$AIC=-101). This significant improvement strongly supported the new version of the model developed in this study and for which we present the estimation results below.

First, we estimate: $\Delta\alpha = 9.31 \times 10^{-3} \text{ h}^{-1}$ and $Q^{VL} = 4.60$ (see Fig. 3 for the latter). Point estimates of the rates of prophage reactivation (expressed in $\text{h}^{-1}$) are thus: $\hat{\alpha}_w = 2.58 \times 10^{-3}$ and $\hat{\alpha}_m = 1.19 \times 10^{-2}$.

Second, we perform non-linear optimizations to estimate the remaining parameters. Overall, computed dynamics of the logit-prevalence and logit-frequencies of the virulent phage closely fit experimental data (Fig. 4). Best MLE estimates are listed in Table 2,

along with their 95% bootstrap-based CIs (see Supplementary Fig. S8 for the density distributions of estimated parameters). We show pairwise correlation coefficients in Supplementary Fig. S9. Some parameters are correlated positively, in particular: $\alpha_w$ with $\alpha_m$, $p_0$ with both reactivation rates, $\phi_m$ with $P_0^{\text{epidemic}}$ and $r$ with $b$; or negatively, in particular: $\tau$ with $b$ and $r$. In Supplementary Fig. S10, we plot the 95% distributions of fitted values. To test the impact of the fixed value used for the burst size B, we study how perturbations of the original value ($80 \text{ virus cell}^{-1}$) affect the estimation of the other parameter values (new point estimates from non-linear optimizations). In Supplementary Fig. S11, we show that all parameters are extremely robust to these perturbations, except for parameter $b$ as expected from Supplementary Fig. S6, (but the product $b \times B$ is always around 3.95). It is also worth noting that the basic reproduction number is not affected by the choice of the value of B, as it always appears as a product with $b$ in (4). We estimate the basic reproduction number of the wild-type strain $\mathcal{R}_{0,w}$ to be 1.48 (95% CI [1.04, 2.15]) and of the virulent strain $\mathcal{R}_{0,m}$ to be 2.20 (95% CI [1.58, 3.05]). Interestingly, the basic reproduction number of the virulent mutant $\mathcal{R}_{0,m}$ is higher than that of the wildtype $\mathcal{R}_{0,w}$, but the virulent mutant always loses the competition in the long term. The basic reproduction number gives the expected number of secondary infections in an otherwise fully susceptible host population (Anderson and May, 1991). The generation interval distribution of the two strains being similar for the lytic cycle, the basic reproduction number provides a good predictor of the relative fitness of the two strains at the early stage of the epidemic (Wallinga and Lipsitch, 2007; Park et al., 2019): the mutant outcompetes the wildtype as $\mathcal{R}_{0,m} > \mathcal{R}_{0,w}$. When the density of susceptible host cells drops, however, the basic reproduction number becomes a poor predictor of fitness. Eventually, only the strain that tolerates the worst environment in terms of resource density—in this case, the lowest density of susceptible host cells—survives (pessimization principle) (Diekmann, 2004). The wildtype replicates better at lower S cells densities than
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Fitted values on experimental data. Fitted values (thick dark lines) were obtained from the simulation based on the best MLE estimations (see Table 2). We estimate model parameters using experimental data (light-colored thin lines) from the evolution experiment designed in Berngruber et al. (2013); initial prevalence is either low (epidemic condition, in red) or high (endemic condition, in blue). Light-colored thick lines correspond to mean logit-values across all replicates per treatment. (a) Logit-prevalence of infection; (b) logit-frequency of cells infected (involved either in a lysogenic or a lytic cycle) by the the mutant (virulent) phage; (c) logit-frequency of the virulent phage in the culture medium (free virus stage).</figcaption></figure>
<table><caption>**Table 2.** MLE parameter estimations. In a first step, parameters $\alpha_w$ and $\alpha_m$ are estimated directly using analytic approximations and linear models. In a second step, we fix $\alpha_w$ and $\alpha_m$ to their previous point estimates and, starting from 2000 sets of initial values, non-linear optimizations maximizing the overall log-likelihood with the experimental data are run to estimate the remaining parameters. 95% CIs are based on Sieve bootstrap (6686 sets of estimates). Throughout, fixed parameters are: $K = 10^9$ cells, $\delta = 0.8 h^{-1}$ and $B = 80$ virus.cell$^{-1}$.</caption><thead><tr><th>Parameter</th><th>Estimated value</th><th>95% CI</th><th>Unit</th></tr></thead><tbody><tr><td colspan="4" style="text-align: center; font-weight: bold;">Parameters of the phage</td></tr><tr><td>α<sub>w</sub></td><td>2.58 × 10<sup>-3</sup></td><td>[1.94 × 10<sup>-4</sup>, 7.87 × 10<sup>-3</sup>]</td><td>h<sup>-1</sup></td></tr><tr><td>α<sub>m</sub></td><td>1.19 × 10<sup>-2</sup></td><td>[1.03 × 10<sup>-3</sup>, 2.85 × 10<sup>-2</sup>]</td><td>h<sup>-1</sup></td></tr><tr><td>φ<sub>w</sub></td><td>0.347</td><td>[0.198, 0.501]</td><td>—</td></tr><tr><td>φ<sub>m</sub></td><td>2.91 × 10<sup>-2</sup></td><td>[7.70 × 10<sup>-3</sup>, 6.94 × 10<sup>-2</sup>]</td><td>—</td></tr><tr><td>a</td><td>1.00 × 10<sup>-6</sup></td><td>[4.34 × 10<sup>-7</sup>, 1.00 × 10<sup>-6</sup>]</td><td>h<sup>-1</sup>.cell<sup>-1</sup></td></tr><tr><td>b</td><td>4.94 × 10<sup>-2</sup></td><td>[2.44 × 10<sup>-2</sup>, 7.51 × 10<sup>-2</sup>]</td><td>—</td></tr><tr><td>τ</td><td>1.08</td><td>[0.83, 3.95]</td><td>h<sup>-1</sup></td></tr><tr><td colspan="4" style="text-align: center; font-weight: bold;">Parameter of the bacteria</td></tr><tr><td>r</td><td>1.60</td><td>[0.77, 5.00]</td><td>h<sup>-1</sup></td></tr><tr><td colspan="4" style="text-align: center; font-weight: bold;">Initial conditions</td></tr><tr><td>p<sub>0</sub><sup>epidemic</sup></td><td>4.24 × 10<sup>-3</sup></td><td>[1.00 × 10<sup>-3</sup>, 1.13 × 10<sup>-2</sup>]</td><td>—</td></tr><tr><td>p<sub>0</sub><sup>endemic</sup></td><td>99.61 × 10<sup>-2</sup></td><td>[0.988, 0.998]</td><td>—</td></tr><tr><td>p<sub>0</sub></td><td>0.416</td><td>[0.400, 0.517]</td><td>—</td></tr><tr><td colspan="4" style="text-align: center; font-weight: bold;">SD of measurement errors</td></tr><tr><td>σ<sub>P</sub></td><td>1.51</td><td>[1.13, 2.00]</td><td>—</td></tr><tr><td>σ<sub>g</sub></td><td>0.52</td><td>[0.52, 0.75]</td><td>—</td></tr><tr><td>σ<sub>q</sub></td><td>0.83</td><td>[0.74, 1.16]</td><td>—</td></tr></tbody></table>

the virulent strain because it relies more on vertical transmission and less to transmission to new susceptible cells. $R_0$ is not maximized in the long-term because the 'niche' of the virus is multidimensional, again because of the two modes of replication (Lion and Metz, 2018).

We compared our estimates of model parameters with previous studies (see reported values in Supplementary Table S3), focusing in particular on the parameters of the phage. Note that previous estimations of these life-history traits were obtained from experimental assays using monomorphic phage cultures. Our estimate of the rate of prophage reactivation of the wildtype strain $\alpha_w$ falls in the range of previous in vitro estimates, which is very large [between $10^{-7}$ and $10^{-2} h^{-1}$ (De Paepe et al., 2016; Little et al., 1999; Zong et al., 2010)], probably due to the variability in experimental methods. For the ci variant, the orders of magnitude of $\alpha_m$ are more consistent (Zong et al., 2010; De Paepe et al., 2016). Other studies, using single-cell monitoring of infected cells to estimate the probability of lysogenization of the wildtype strain $\phi_w$, obtained very similar results when the multiplicity of infection (MOI) is low, around 0.35 on average (Zeng et al., 2010)—but see the section 'Discussion' for the effect of MOI. We estimate a 12-fold decrease for the probability of lysogenization of the virulent strain $\phi_m$ and, to our knowledge, this parameter has not been estimated elsewhere. The lysis time $1/\tau$, estimated close to 1 h (though with larger uncertainty), is also consistent with the biology of phage $\lambda$ (Zong et al., 2010; De Paepe and Taddei, 2006; De Paepe et al., 2016; Lindberg et al., 2014; Shao and Wang, 2008). On the other hand, the adsorption rate $a$ is poorly estimated. Although this issue did not arise with simulated data, our estimates based on experimental data are mostly stuck at the upper bound ($10^{-6}$), several orders of magnitude higher than expected (De Paepe and Taddei, 2006; De Paepe et al., 2016; Lindberg et al., 2014; Shao and Wang, 2008).

Wang, 2008). To investigate the sensitivity of the other estimates to the estimated value of $a$, we reiterate non-linear optimizations fixing $a$ to different values ranging from $10^{-9}$ to 1 to compute new point estimates. The inference of most parameters was robust to the value of $a$, with no relative variation higher than 50% found for values of $a$ between $10^{-8}$ to 1 (Supplementary Fig. S12A). We carry out the same sensitivity analysis for the intrinsic growth rate of the host $r$ which was also poorly estimated (its distribution being mainly flat): with values of $r$ ranging between 1 and 5, the probabilities of lysogenization are little affected by the perturbations in $r$ (relative variation contained between -7% and +15%) while parameters $b$ and $\tau$—the most correlated with $r$—are the most sensitive (Supplementary Fig. S12B).

## 4. Discussion

In a previous study, Berngruber et al. (2013) carried out an evolution experiment with the phage $\lambda$ to validate several theoretical predictions on the evolution of virulence and transmission. This study used a two-step approach: (i) a mathematical model tailored to the biology of phage $\lambda$ was used to show how epidemiological dynamics feedback on the evolution of the pathogen and (ii) tracking the variation of the frequency of a viral mutant across time and across compartments (infected host and free virus) confirmed the impact of the density of susceptible cells on the transient evolution of the virus at different stages of its life cycle. This experimental validation of the theory focused on the qualitative match of the experimental results with the predictions of the model and provided empirical support for evolutionary epidemiology theory. In the present study, we adopt a reverse approach where we start from the data and try to improve the theoretical model developed to describe the evolutionary epidemiology of phage $\lambda$. We calculate the selection gradients and the differentiation across the life stages of the virus at the onset and at the end of the epidemic. This analysis allows us to better quantify the processes driving the evolutionary dynamics of the virus population.

First, we noticed that the previous model failed to properly capture the transient evolutionary dynamics in the infected compartment. We improved the goodness of fit with the experimental data by distinguishing two types of infected cells: lysogens (L) and cells committed to the lytic pathway (Y). Indeed, adding extra stages allowed us to observe a transitory peak in the frequency of hosts infected by the virulent phage (Supplementary Fig. S13A, middle panel), similar to what we observed in the experimental data. In this model, the peak is due to the frequency of the virulent phage among Y bacteria that, unlike L bacteria, transiently overshoots during the acute phase of the epidemic. Crucially, lysis was assumed to be instantaneous in the original model (Berngruber et al., 2013) and including this additional stage in the phage life cycle allowed us to take the lag between infection and lysis into account (Yin and McCaskill, 1992; You and Yin, 1999; Mitarai et al., 2016; Brown et al., 2022; Geng et al., 2023; De Paepe et al., 2016). For the sake of simplicity and parsimony, we only kept a single extra stage Y, yielding exponentially distributed lysis time (Mitarai et al., 2016).

Second, we performed a theoretical analysis of this new model to go beyond the numerical approach used in Berngruber et al. (2013). In particular, we obtained analytic approximations for the spread of the virulent phage at the beginning and at the end of the epidemic: (i) the virulent mutant increases in frequency at the onset of the epidemic (when the density of susceptible hosts is high) with a speed approximately proportional to $-\Delta\phi$; (ii) the virulent mutant decreases in frequency at the end of the epidemic
(when the density of susceptible hosts is low) with a negative selection gradient –$\Delta\alpha$, which is consistent with the prediction that in the long-term, and in the absence of an influx of susceptible hosts, the temperate phage should evolve a fully lysogenic strategy where $\alpha=0$ and $\phi=1$ (Bruce et al., 2021; Wahl et al., 2019); and (iii) the virulent phage is always more frequent in the free virus stage than among prophages and this differentiation is approximated by $1+\Delta\alpha/\alpha_w$ at the end of the epidemic. The first approximation only held for a very limited period of time because the density of susceptible cells drops very rapidly and this epidemiological feedback affects the selection on transmission and virulence. In contrast, the final two predictions were valid as soon as the epidemic reaches high prevalence and, interestingly, yielded a novel way to estimate the rates of prophage reactivation of both strains.

Third, we developed a statistical inference approach to estimate the parameters of the model. We implemented a two-step MLE procedure: (i) we first obtained point estimates of the reactivation rates $\alpha_w$ and $\alpha_m$ of both viral strains and (ii) we then fixed $\alpha_w$ and $\alpha_m$ to their point estimates and ran non-linear optimizations to infer the remaining parameters of the model (except the burst size $B$ which had to be fixed because of an identifiability issue). The Sieve bootstrap technique was used to generate joint distributions for all the estimated parameters (Table 2). We showed that our estimates of the key phenotypic traits $\alpha_w$, $\alpha_m$, and $\phi_w$ were consistent with existing literature (Supplementary Table S3).

Crucially, our inference approach is quite different from previous studies. First, we analyse the dynamics of a polymorphic viral population. Second, we use three types of data to infer the parameter values of the model: (i) epidemiological data (i.e. the prevalence of the infection), (ii) temporal variations in variant frequencies, and (iii) differentiation of the variant frequency across compartments. Each data type carries complementary information; together, they allow us to jointly estimate the life-history traits of both strains of the phage. This novel method is particularly well suited to estimate the rates of prophage reactivation, for which only the endemic treatment is needed—see Supplementary Appendix S3.2 where we propose a Bayesian counterpart to easily obtain posterior distributions and credible intervals for these two parameters (applied to our experimental data, we show results in Supplementary Fig. S14A). Estimating the probabilities of lysogenization is however more difficult as it requires to monitor the transient dynamics of the epidemic treatment. The rapid epidemiological feedback that occurred in emerging epidemics makes it necessary to consider both the epidemiology and the evolution of the infection. It might thus be relevant to carry out shorter experiments for the epidemic treatment—i.e. focusing on the early state of the epidemic—but with a more frequent sampling effort to monitor the change in frequency more precisely. Alternatively, the use of two-stage chemostats, where the influx of susceptible hosts could be maintained experimentally (Bonachela and Levin, 2014; Husimi et al., 1982), might also be an option worth investigating.

Even though we have improved the original model, many features of the present model could be challenged. For example, we assume all phenotypic traits to be constant across time and across chemostats. Yet, key life-history traits of phage $\lambda$ are known to vary with the environment (Wegrzyn and Wegrzyn 2005) (i.e. phenotypic plasticity). In particular, the probability of lysogenization and the reactivation rate of the $\lambda$ mutant used in the evolution experiment is temperature-sensitive (Berngruber et al., 2013; St-Pierre and Endy, 2008; Zong et al., 2010). Interestingly, the estimated rates of reactivation varied among chemostats (Supplementary Fig. S14B, using a Bayesian approach). This variation

may have been driven by small variations in temperature among the chemostats that could have affected the reactivation rate of the mutant (higher rates of reactivation would be expected with higher temperature). The probability of lysogenization of phage $\lambda$ is also known to be more likely to occur at high MOI (Kourilsky, 1973; Kourilsky and Knapp, 1975; Zeng et al., 2010). MOI-dependent phenotypic plasticity may also affect the rate of reactivation of temperate phages (Bruce et al., 2021), or the lysis time of virulent phages (Rutberg and Rutberg, 1965). Yet, our model has the advantage of the parsimony, while already convincingly reproducing the qualitative patterns observed in the data. More sophisticated models would improve the fit to the data compared to simpler, nested models. Additional details may be rooted in experimental knowledge on the system, but the data may not be rich enough to infer the extra-parameters precisely (and accordingly, these more complex models may not be selected over simpler ones in statistical model comparisons). Moreover, simpler models are easier to analyse mathematically and lend themselves more easily to interpretation. It is thus delicate to know where to draw the line between the parsimony and the goodness of fit of a model. Navigating in this trade-off, a first model was developed capturing the transient increase in virulent strain (Berngruber et al., 2013), and we extended it with the known delay between phage infection and lysis to better capture the similarity in variant frequency in cells and free viruses. Our final model fits the data well, but not all parameters can be inferred precisely, in ways that depend on the details of the model and would have been difficult to anticipate before formally fitting the model.

Our study emphasizes the benefits of combining theoretical, statistical, and experimental approaches. Combining these perspectives effectively requires an iterative process. Each step is an opportunity to challenge and improve our understanding of the biological model, to elucidate microbial life cycles and to provide support and guidance during experimentation. While most experimental life-history assays only focus on monomorphic population, we demonstrate the relevance of evolution experiments where different strains are put in competition and where tracking their relative frequencies over time may yield novel and efficient ways to estimate some dimensions of their phenotypes, especially when it is difficult to access absolute densities. This approach could be used under different conditions and in particular *in vivo* (De Paepe et al., 2016). The present study may thus help characterize the phenotypic traits of microbial organisms used in experimental evolution. More broadly, this overall quantitative process is similar in spirit with what can be done in epidemiological studies, particularly in public health—though generally with a much simpler description of the life cycle. During the COVID-19 pandemic, for instance, the successive emergence and sweeps of SARS-CoV-2 variants of concern has raised many questions about the short- and long-term evolution of the virus, especially in terms of transmission and virulence. Statistical inference based on demographic/epidemiological data (e.g. prevalence, deaths, hospital admissions) and genetic data (e.g. strain frequencies) quantified variants properties such as transmission, virulence, infectious period, immune escape, or generation interval distribution [e.g. Davies et al. (2021); Blanquart et al. (2022); Benhamou et al. (2023)]. In many scenarios, the frequency of a variant could be measured in different compartments, for example between naive and vaccinated hosts (e.g. for SARS-CoV-2 variants Omicron), or between infected hosts and an environmental compartment (e.g. enteropathogenic E. coli in environmental waters). As shown in this study, interpreting the temporal variation in the selection gradient and the differentiation of a novel
variant between distinct compartments may be a powerful tool to clarify the relationship between its frequencies in distinct life stages and to better characterize its phenotype.

## Supplementary data

Supplementary data are available at Virus Evolution online.

**Conflict of interest:** None declared.

## Funding

We acknowledge the French Ministry of Higher Education and Research and the École Normale Supérieure Paris-Saclay for the PhD scholarship of WB.

## Data availability

Data that were used in this study, along with the scripts for the analyses (in R), are available in this GitHub repository as well as at https://zenodo.org/records/13138349.

## References

* Alizon S, Hurford A, Mideo N, et al. Virulence evolution and the trade-off hypothesis: history, current state of affairs and the future. *J Evol Biol* 2009; **22**:245–259.

* Alizon S, Michalakis Y. Adaptive virulence evolution: the good old fitness-based approach. *Trends Ecol Evol* 2015; **30**:248–254.

* Anderson RM, May RM. Coevolution of hosts and parasites. *Parasitology* 1982; **85**:411–426.

* Anderson RM, May RM. *Infectious Diseases of Humans: Dynamics and Control*. Oxford University Press, 1991.

* Benhamou W, Lion S, Choquet R, et al. Phenotypic evolution of SARS-CoV-2: a statistical inference approach. *Evolution* 2023; **77**:2213–2223.

* Berngruber TW, Froissart R, Choisy M, et al. Evolution of virulence in emerging epidemics. *PLoS Pathogens* 2013; **9**:e1003209.

* Berngruber TW, Weissing FJ, Gandon S. Inhibition of superinfection and the evolution of viral latency. *J Virol* 2010; **84**:10200–10208.

* Blanquart F, Hozé N, Cowling BJ *et al*. Selection for infectivity profiles in slow and fast epidemics, and the rise of SARS-CoV-2 variants. *Elife* 2022; **11**:e75791.

* Bonachela JA, Levin SA. Evolutionary comparison between viral lysis rate and latent period. *J Theor Biol* 2014; **345**:32–42.

* Brown S, Mitarai N, Sneppen K. Protection of bacteriophage-sensitive *Escherichia coli* by lysogens. *Proc Natl Acad Sci*, 2022, **119**:e2106005119.

* Bruce JB, Lion S, Buckling A *et al*. Regulation of prophage induction and lysogenization by phage communication systems. *Curr Biol* 2021; **31**:5046–5051.

* Bühlmann P. Sieve bootstrap for time series. *Bernoulli* 1997; 123–148.

* Bull JJ. Virulence. *Evolution* 1994; **48**:1423–1437.

* Davies NG, Abbott S, Barnard RC, *et al*. Estimated transmissibility and impact of SARS-CoV-2 lineage B. 1.1.7 in England. *Science* 2021; **372**:eabg3055.

* Day T. On the evolution of virulence and the relationship between various measures of mortality. *Proceedings of the Royal Society of London. Series B: Biological Sciences*, **269**, 1317–1323, 2002.

* Day T, Proulx SR. A general theory for the evolutionary dynamics of virulence. *Am Naturalist* 2004; **163**:E40–E63.

* De Paepe M, Taddei F. Viruses' life history: towards a mechanistic basis of a trade-off between survival and reproduction among phages. *PLoS Biol* 2006; **4**:e193.

* De Paepe M, Tournier L, Moncaut E *et al*. Carriage of λ latent virus is costly for its bacterial host due to frequent reactivation in monoxenic mouse intestine. *PLoS Genet* 2016; **12**:e1005861.

* Diekmann O. *A Beginner's Guide to Adaptive Dynamics*, Vol. 63, Banach Center Publications, 2004, 47–86.

* Diekmann O, Heesterbeek J, Roberts MG. The construction of next-generation matrices for compartmental epidemic models. *J R Soc Interface* 2010; **7**:873–885.

* Echols H. Developmental pathways for the temperate phage: lysis vs lysogeny. *Annu Rev Genet* 1972; **6**:157–190.

* Frank SA. Models of parasite virulence. *Q Rev Biol* 1996; **71**:37–78.

* Gandon S. Why be temperate: lessons from bacteriophage λ. *Trends Microbiol* 2016; **24**:356–365.

* Gandon S, Day T. The evolutionary epidemiology of vaccination. *J R Soc Interface* 2007; **4**:803–817.

* Geng Y, Nguyen TVP, Homaee E, *et al*. Using population dynamics to count bacteriophages and their lysogens. *bioRxiv* 2023; 2023–10.

* Husimi Y, Nishigaki K, Kinoshita Y, et al. Cellstat—a continuous culture system of a bacteriophage for the study of the mutation rate and the selection process at the DNA level. *Rev Sci Instrum* 1982; **53**:517–522.

* Hyndman RJ, Khandakar Y. Automatic time series forecasting: the forecast package for r. *J Stat Softw* 2008; **27**:1–22.

* Kelley CT. *Iterative Methods for Optimization*. SIAM, 1999.

* Kourilsky P. Lysogenization by bacteriophage lambda: I. multiple infection and the lysogenic response. *Mol Gen Genet* 1973; **122**:183–195.

* Kourilsky P, Knapp A. Lysogenization by bacteriophage lambda: III. Multiplicity dependent phenomena occurring upon infection by lambda. *Biochimie* 1975; **56**:1517–1523.

* Lenski RE, May RM. The evolution of virulence in parasites and pathogens: reconciliation between two competing hypotheses. *J Theor Biol* 1994; **169**:253–265.

* Lindberg HM, McKean KA, Wang I-N. Phage fitness may help predict phage therapy efficacy. *Bacteriophage* 2014; **4**:e964081.

* Lion S, Metz JA. Beyond R0 maximisation: on pathogen evolution and environmental dimensions. *Trends Ecol Evol* 2018; **33**:458–473.

* Little JW, Shepley DP, Wert DW. Robustness of a gene regulatory circuit. *EMBO J* 1999; **18**:4299–4307.

* Lwoff A. Lysogeny. *Bacteriol Rev* 1953; **17**:269–337.

* Mitarai N, Brown S, Sneppen K. Population dynamics of phage and bacteria in spatially structured habitats using phage λ and *Escherichia coli*. *J Bacteriol* 2016; **198**:1783–1793.

* Nelder JA, Mead R. A simplex method for function minimization. *Comput J* 1965; **7**:308–313.

* Park SW, Champredon D, Weitz JS, et al. A practical generation-interval-based approach to inferring the strength of epidemics from their speed. *Epidemics* 2019; **27**:12–18.

* Ptashne M. A Genetic Switch: Phage λ and Higher Organisms. Oxford: Blackwell Publishers, 1992.

* R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing, 2022.

* Rutberg B, Rutberg L. Role of superinfecting phage in lysis inhibition with phage t4 in *Escherichia coli*. *J Bacteriol* 1965; **90**:891–894.

* Shao Y, Wang I-N. Bacteriophage adsorption rate and optimal lysis time. *Genetics* 2008; **180**:471–482.

* Soetaert K, Petzoldt T, Setzer RW. Solving differential equations in R: package deSolve. *J Stat Softw* 2010; **33**:1–25.

* St-Pierre F, Endy D. Determination of cell fate selection during phage lambda infection. *Proc Natl Acad Sci*, 2008; **105**:20705–20710.

* Sussman R, Jacob F. On a thermosensitive repression system in the *Escherichia coli* lambda bacteriophage. *Comptes Rendus*
héBdomadaires des séAnces de l'Académie des sciences 1962; **254**:
1517–1519.

Ulloa G, Allende-Cid H, Allende H. Sieve bootstrap prediction intervals for contaminated non-linear processes. In *Progress in Pattern Recognition, Image Analysis, Computer Vision, and Applications: 18th Iberoamerican Congress, CIARP 2013, Proceedings, Part I*, **18**, pp.84–91.
Springer: Havana, Cuba, 2013.

Wahl LM, Betti MI., Dick DW et al. Evolutionary stability of the lysis-
lysogeny decision: why be virulent?. *Evolution* 2019; **73**:92–98.

Wallinga J, Lipsitch M. How generation intervals shape the relation-
ship between growth rates and reproductive numbers. *Proc R Soc
B: Biol Sci*, 2007; **274**:599–604.

Wang I-N. Lysis timing and bacteriophage fitness. *Genetics* 2006; **172**:17–26.

Wegrzyn G, Wegrzyn A. Genetic switches during bacteriophage λ development. *Prog Nucleic Acid Res Mol Biol* 2005; **79**: 1–48.

Yin J, McCaskill J. Replication of viruses in a growing plaque: a reaction-diffusion model. *Biophys J* 1992; **61**: 1540–1549.

You L, Yin J. Amplification and spread of viruses in a growing plaque.
*J Theor Biol* 1999; **200**:365–373.

Zeng L, Skinner SO, Zong C *et al*. Decision making at a subcellular level determines the outcome of bacteriophage infection. *Cell* 2010; **141**:682–691.

Zong C, So L, Sepúlveda LA et al. Lysogen stability is determined by the frequency of activity bursts from the fate-determining gene. *Mol Syst Biol* 2010; **6**:440.