Assessing the impact of bacterial heterogeneity on bacteriophage population dynamics

TL;DR

This study examines how variability among bacteria affects how bacteriophages grow and spread, showing that traditional models may not capture this complexity accurately.

Contribution

The paper introduces a new modeling framework that accounts for cell-level variability in phage-host dynamics.

Findings

Only the distribution of the latent period and its effect on burst size significantly impacts population dynamics.

Classical models work well at low variability but fail to capture dynamics at high variability.

The new distributed framework is more accurate when high variability is present.

Abstract

Bacteriophages follow complex replication dynamics in codependence with their host cells. The corresponding dynamics are usually described with three key parameters: the adsorption rate \documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\varphi$$\end{document}, the burst size \documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\beta$$\end{document} and the latent period \documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$\tau$$\end{document}. Single-cell experiments have shown cellular variability in latent period and burst size. Yet, those parameters are traditionally modelled as uniform values for the whole bacterial population, which may introduce bias in the predicted profiles. Here, we systematically assessed the relevance of modelling this cell-level variability in the context of bacteriophage dynamics. We developed a comprehensive modelling framework, aiming to: (i) identify which parameter distributions impact population dynamics, (ii) quantify how difference in distribution can influence predictions, and (iii) inform model selection by comparing relevant modelling approaches. A distributed delay differential equation model with randomly distributed latent period, burst size and adsorption rate within the cell population was developed. It also included logistic growth of the host cell, saturable adsorption rate and intracellular replication dynamics. This model was compared to classical models, both fixed delay differential equation and transit compartments models, over typical ranges of phage parameters. Only the distribution of the latent period and its influence on burst size impacted population dynamics. The two classical models proved to be good approximations of dynamics, but only at low cellular variability. As such, although the classical approaches can be preferred at lower variability, this distributed framework is warranted when high variability is expected, or when the other approaches fail to accurately capture the dynamics.

The online version contains supplementary material available at 10.1007/s10928-026-10025-y.

Figures

Single-step growth curves simulated with a fixed latent period (red curve) and a distributed latent period (blue curve)

Workflow of the study: A. Definition of parameter distributions to compare (Section “Considered distributions”). B. Simulation of the corresponding bacteria and phage time profiles with a distributed delay differential equation (DDDE) model (Section “Distributed parameter framework”). C. Comparison of resulting bacterial and viral load profiles for different distributions (Section “Evaluation of the impact of the distribution”). D. Compilation of the differences observed after repeating steps A to C for 100 sets of parameters, to quantitatively evaluate the overall effect of parameter distribu

Compartmental representation of the distributed delay differential equations model with intracellular replication dynamics

Single-step growth curves simulated with a fixed latent period (red curve) and a distributed latent period (blue curve)

Workflow of the study: A. Definition of parameter distributions to compare (Section “Considered distributions”). B. Simulation of the corresponding bacteria and phage time profiles with a distributed delay differential equation (DDDE) model (Section “Distributed parameter framework”). C. Comparison of resulting bacterial and viral load profiles for different distributions (Section “Evaluation of the impact of the distribution”). D. Compilation of the differences observed after repeating steps A to C for 100 sets of parameters, to quantitatively evaluate the overall effect of parameter distribu

Compartmental representation of the distributed delay differential equations model with intracellular replication dynamics

Compartmental representation of the transit compartment model, for an Erlang-distributed latent period

Differences in bacterial and viral load predictions between log-normally, normally and gamma distributed latent period for equal distribution width, compiled over 500 random sets of parameters

Error on bacterial and viral load predictions induced by fixing the latent period to early lysis time, to the median latent period, or to the mean latent period, compiled over 500 random sets of parameters. The reference used is a log-normally distributed latent period

Error on bacterial and viral load predictions of the DDE and TC models relative to the full DDDE model, compiled over 500 random sets of parameters. The reference used is a log-normally distributed latent period

Full text

Introduction

Bacteriophages or phages have gained renewed interest in the past decade due to the increasing threat of multidrug-resistant bacterial infections and their potential utility as anti-infective therapy [1]. Bacteriophages are viruses which infect target bacteria and exhibit complex population dynamics in interaction with their hosts. In order to optimise bacteriophage therapy, we need to quantitatively understand and predict the effects of said interactions. To that end, mathematical models describing bacteria-phage interactions have been previously developed [2, 3]. These models typically consider three key parameters: (i) the rate at which phages bind to the target bacteria (adsorption rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi$$\end{document} ); (ii) the delay between infection and lysis of the host bacterium (latent period \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} ); and (iii) the number of new phages released per lysed bacterium (burst size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} ). While such models provide valuable insights, they generally neglect the intrinsic heterogeneity between individual bacteria, assuming uniform parameter values across the entire bacterial population [2].

Single-cell studies have demonstrated substantial variability in both latent period [4] and burst size [5] between individual bacteria. Moreover, experimental studies have shown that even under conditions of synchronised infection [6], phage concentrations increase gradually over time rather than abruptly [7], which is consistent with a distributed latent period (Fig. 1). Yet, the latent period is typically reported as a uniform value, when the first lysis is observed experimentally [8]. This approach not only disregards the parameter distribution, but also could lead to a systematic bias with regard to the latent period observed at the single-cell level [9].Fig. 1. Single-step growth curves simulated with a fixed latent period (red curve) and a distributed latent period (blue curve)

The distributed nature of the latent period has been described using distributed delay differential equations (DDDE) [10, 11]. However, such models are rarely used in practice. Their limited adoption is largely due to their theoretical complexity, long computation times, and the scarcity of accurate experimental data on relevant parameter distributions. It remains unclear whether and under which conditions incorporating cellular variability meaningfully improves model predictions of population dynamics.

The current study aimed to address this gap by: i) identifying analytically parameter distributions which impact population dynamics; ii) evaluating the sensitivity of population dynamics to the shapes of these distributions; (iii) quantifying the impact of fixing the parameters to representative values on population dynamics; and (iv) informing model selection by comparing relevant modelling alternatives. To this end, we used a theoretical bottom-up modelling approach. We first developed a bacteriophage pharmacodynamic model accounting for cellular variability (Sections “Base model structure” & “Distributed parameter framework”). We then used this model to simulate bacteria and phage dynamics (Fig. 2B) for different parameter distributions (Section “Considered distributions”, Fig. 2A). We compared the obtained profiles and defined an error metric to quantify their differences (Fig. 2C). To systematically and quantitatively assess how the parameter distributions impact the dynamics, we repeated the previous steps with 100 different sets of model parameters (Section “Evaluation of the impact of the distribution”, Fig. 2D).Fig. 2. Workflow of the study: A. Definition of parameter distributions to compare (Section “Considered distributions”). B. Simulation of the corresponding bacteria and phage time profiles with a distributed delay differential equation (DDDE) model (Section “Distributed parameter framework”). C. Comparison of resulting bacterial and viral load profiles for different distributions (Section “Evaluation of the impact of the distribution”). D. Compilation of the differences observed after repeating steps A to C for 100 sets of parameters, to quantitatively evaluate the overall effect of parameter distribution on population dynamics (Section “Evaluation of the impact of the distribution”)

Theoretical

Base model structure

The study was based on a classic mathematical model for phage dynamics [2]:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{dU}{dt}(t)&= r(t) \cdot U(t) - \delta _{U} \cdot U(t) \nonumber \\&\qquad \qquad \,\,\,- \varphi (t) \cdot U(t) \cdot P(t) \nonumber \\ \frac{dI}{dt}(t)&= \varphi (t) \cdot U(t) \cdot P(t) - \delta _I \cdot I(t) \\&\qquad \qquad \,\,\,- e^{- \delta _I \tau } \cdot \varphi (t-\tau ) \cdot U(t-\tau ) \cdot P(t-\tau ) \nonumber \\ \frac{dP}{dt}(t)&= \beta \cdot e^{- \delta _I \tau } \cdot \varphi (t-\tau ) \cdot U(t-\tau ) \cdot P(t-\tau ) \nonumber \\&\qquad \qquad \,\,\,- \delta _P \cdot P(t) - \varphi (t) \cdot (U(t) + I(t)) \cdot P(t)\nonumber \end{aligned}$$\end{document}

where U, I and P are the concentrations of uninfected bacteria, infected bacteria, and phages respectively, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _U, \delta _I, \delta _P$$\end{document} are their respective decay rates, r(t) is the bacterial growth rate, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi (t)$$\end{document} is the adsorption rate, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} is the burst size and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} is the latent period. Here, phages may adsorb to either uninfected or infected bacteria with the adsorption rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi (t)$$\end{document} .

Equation 1 depends on the states of the variables at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t-\tau$$\end{document} , which makes this a delay differential equations (DDE) model. Such equations require information on the history of the state variables. Here, t=0 is taken as the time of phage administration, hence:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \forall t < 0, P(t) = 0 \end{aligned}$$\end{document}

Any history of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi (t)$$\end{document} and U(t) thus results in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\forall t < 0, \varphi (t) \cdot P(t) \cdot U(t) = 0$$\end{document} . This history of P, as well as usual initial conditions at t=0 for the state variables, is consequently enough to fully determine the system.

The bacterial growth was assumed logistic [12], which gives:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} r(t) = r_{max} \cdot (1 - \frac{U(t) + I(t)}{K_C}) \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{max}$$\end{document} is the max bacterial growth rate and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_{C}$$\end{document} is the ceiling concentration of bacterial growth.

A saturable adsorption rate was assumed as follows [13]:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varphi (t) = \varphi _{max} \cdot \frac{P_{50}}{P_{50} + P(t)} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _{max}$$\end{document} is the maximal adsorption rate and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_{50}$$\end{document} is the phage concentration at which the adsorption rate is half its maximal value.

Distributed parameter framework

To account for cellular variability, Eq. 1 was adapted with distributed parameters. Prior cell-studies showed evidence of dependence of burst size on latent period [14], which was accounted for by introducing conditional distributions.

The latent period \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} is randomly distributed around a fixed mean value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _0$$\end{document} . The corresponding probability density is denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{\tau }(\tau )$$\end{document} . A minimal possible latent period \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{min}$$\end{document} can also be defined to account for the fact that a phage cannot instantly lyse a bacterium [10].

For a given value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} , the burst size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} is randomly distributed around a mean value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _0(\tau )$$\end{document} , for which the conditional probability density is written as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{\beta |\tau }(\beta | \tau )$$\end{document} . The relationship between mean burst size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _0(\tau )$$\end{document} and latent period \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} was described using the following phenomenological expression [14]:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \beta _0(\tau ) = \beta _{max} \frac{e^{r_\beta (\tau - D)} - 1}{e^{r_\beta (\beta \tau _{50} - D)} + e^{r_\beta (\tau - D)} - 2} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{max}$$\end{document} is the maximal burst size, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{\beta }$$\end{document} the intracellular phage replication rate, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta \tau _{50}$$\end{document} is the latent period corresponding to half of the maximum burst size and D is a time delay below which burst size is 0. A fixed value for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _0(\tau )$$\end{document} , independent of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} , was also explored (appendix).

The population mean burst size, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{0,mean} = \int _0^{+\infty } \beta _0(\tau )f_\tau (\tau )d\tau$$\end{document} , is the parameter typically determined from single-step growth curves [8]. Therefore, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{max}$$\end{document} in Eq. 5 was re-expressed as a function of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{0,mean}$$\end{document} (derivation in appendix):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \beta _{max} = \frac{\beta _{0, mean}}{\int _{\tau = \tau _{min}}^{+\infty } \frac{e^{r_\beta (\tau - D)} - 1}{e^{r_\beta (\beta \tau _{50} - D)} + e^{r_\beta (\tau - D)} - 2} \cdot f_\tau (\tau ) \cdot d\tau } \end{aligned}$$\end{document}

The adsorption rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi (t)$$\end{document} is randomly distributed around a mean value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _{0}(t)$$\end{document} that is a function of time at which the phage adsorbs to the bacteria (assumed to follow Eq. 4). The corresponding conditional probability density is written \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{\varphi |t}(\varphi |t)$$\end{document} .

The associated differential equations are:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{dU}{dt}(t)&= r(t) \cdot U(t) \nonumber \\&\qquad \quad - \int _{\varphi = 0}^{+ \infty } \varphi \cdot U(t) \cdot P(t) \cdot f_{\varphi | t}(\varphi | t) \cdot d\varphi - \delta _{U} \cdot U(t)\nonumber \\ \frac{dI}{dt}(t) =&\int _{\varphi = 0}^{+ \infty } \varphi \cdot U(t) \cdot P(t) \cdot f_{\varphi | t}(\varphi | t) \cdot d\varphi \\&- \int _{\tau = 0}^{t} \int _{\varphi = 0}^{+ \infty } e^{- \delta _I \tau } \cdot \varphi \cdot U(t-\tau ) \cdot P(t-\tau ) \nonumber \\&\cdot f_{\varphi | t}(\varphi | t - \tau ) \cdot f_{\tau }(\tau ) \cdot d\varphi \cdot d\tau - \delta _I \cdot I(t)\nonumber \\ \frac{dP}{dt}(t) =&\int _{\tau = 0}^{t} \int _{\varphi = 0}^{+ \infty } \int _{\beta = 0}^{+ \infty } \beta \cdot e^{- \delta _I \tau } \cdot \varphi \cdot U(t-\tau ) \cdot P(t-\tau ) \nonumber \\&\qquad \quad \cdot f_{\beta | \tau }(\beta | \tau ) \cdot f_{\varphi | t}(\varphi | t-\tau ) \cdot f_{\tau }(\tau ) \cdot d\beta \cdot d\varphi \cdot d\tau \nonumber \\&- \int _{\varphi = 0}^{+ \infty } \varphi \cdot (U(t) + I(t)) \cdot P(t) \cdot f_{\varphi | t}(\varphi | t) d\varphi - \delta _P \cdot P(t)\nonumber \end{aligned}$$\end{document}

This model (represented in Fig. 3) is a distributed delay differential equation (DDDE) model accounting for the dependence of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} (Eq. 5), associated with intracellular replication dynamics. Those equations are a naive representation of what is referred to as the "Full DDDE" model. Similarly to the DDE model, Eq. 2 provides sufficient history for this model structure. In particular, it is not necessary to integrate for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau> t$$\end{document} because \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\forall \tau> t, P(t-\tau ) = 0.$$\end{document} Fig. 3. Compartmental representation of the distributed delay differential equations model with intracellular replication dynamics

Methods

Evaluation of the impact of the distribution

Analytical assessment

An analytical analysis of the differential (Eq. 7) was first conducted to gain general insights on whether the distributions of the latent period \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} , the burst size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} and the adsorption rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi$$\end{document} influence population dynamics.

The roles of the distributions of the burst size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} and the adsorption rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi$$\end{document} in population dynamics were fully analytically determined. However, the analytical study was insufficient for the latent period \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} . Further simulation work was thus required. As such, the impact of the latent period distribution was investigated through both visual and quantitative analyses of simulated scenarios.

Exploratory analysis

Time profiles of bacteria and bacteriophages were simulated with different latent period distributions. Those profiles were used to inspect how the selection of the distribution model influenced the shape of the simulated curves. In the exploratory study, a visual of the profiles was repeated for various sets of parameters. An example of profiles was included in the results section.

Quantitative analysis

The difference between two simulated profiles A and B was quantified with a metric based on the normalised area between the curves (Fig. 2C):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {Error} (\%) = 100 \cdot \frac{2 \cdot \text {Difference Area}}{AUC_A + AUC_B} = 100 \cdot \frac{2 \cdot \int _{t=0}^{+\infty } |X_A - X_B| \textit{ dt}}{\int _{t=0}^{+\infty } X_A \textit{ dt} + \int _{t=0}^{+\infty } X_B \textit{ dt}} \end{aligned}$$\end{document}

Given any two distribution models, a Monte Carlo analysis of the difference between the two models for both bacteria and phage profiles was conducted. 500 sets of parameters were randomly generated, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{max}, \delta _U, \delta _I, \delta _P, \tau _0, \tau _{min}, r_{\beta }, \varsigma _{eff} \text { and } \beta \tau _{50}$$\end{document} uniformly sampled, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_C, \beta _{max}, \varphi _0, P_{50}, U_{0} \text { and } P_{0}$$\end{document} log-uniformly sampled (to cover their different orders of magnitude). The results were compiled in bar diagrams for both bacterial and viral load predictions.

Considered distributions

Shapes of distribution considered

Different shapes of latent period distribution were simulated to investigate whether the population dynamics were sensitive to it:

• The truncated normal distribution.
• The log-normal distribution.
• The gamma distribution, often used in distributed lifespan models because of its direct equivalence to a transit compartments model when the shape parameter is a positive integer [15]. In each scenario, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$sd_\tau$$\end{document} was defined as the standard deviation of the latent period distribution, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$rsd_\tau = \frac{sd_\tau }{\tau _0}$$\end{document} the associated relative standard deviation. Distributions with relative standard deviations ranging from 10% to 50% were considered to match ranges of variability observed experimentally [4, 16, 17].

Scenarios with fixed parameters

Different scenarios with a non-distributed latent period were also explored, to properly quantify the error associated with neglecting the distribution with the classical DDE model:

• Fixing the latent period to the mean of the associated distribution - referred to as "Fixed Mean".
• Fixing the latent period to the median of the associated distribution - referred to as "Fixed Median".
• Fixing the latent period to time of first observed lysis - referred to as "Fixed Early". The time of first observed lysis was computed as the time after which the number of infective centers has doubled in a single-step growth experiment (see appendix).

Comparison with transit compartments

Transit compartments (TC) models (Fig. 4) provide a computationally efficient method to capture stochasticity in the latent period [9], equivalent to a DDDE model with Erlang-distributed latent period [15]. However, TC models do not capture the potential relationship between burst size and latent period, as presented in Eq. 5, and assume an unvarying mean burst size: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\forall \tau , \beta _0(\tau ) = \beta _{0,mean}$$\end{document} . Moreover, TC models do not account for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{min} \ne 0$$\end{document} . These represent constitutive differences with the full DDDE model, and the two are not strictly equivalent. Thus, the TC model can be seen as an approximation of the full DDDE model, and its performance was compared to the approximation with the DDE model. The number of transit compartments N and the transit rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_{tr}$$\end{document} were defined as [15]:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} N = round(\frac{1}{rsd_\tau ^2}) , k_{tr} = \frac{N}{\tau _0} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$rsd_{\tau }$$\end{document} are respectively the mean and the relative standard deviation of the distribution of the latent period.Fig. 4. Compartmental representation of the transit compartment model, for an Erlang-distributed latent period

Simulation parameters

In the simulations, parameter ranges proximate to typical values reported in literature were used (Table 1), for different bacteria and phage inoculums ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U_{init}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_{init}$$\end{document} respectively).Table 1. Parameter ranges used in the simulationsParameterSimulated rangeReferencesBacterial growth \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{max} \text (/h)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.5 - 2.5$$\end{document} [18] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_{C} \text (bact/mL)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^8 - 10^{10}$$\end{document} [19, 20]Phage adsorption \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _{0} \text (mL / cell / h)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-9} - 10^{-7}$$\end{document} [21–23] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_{50} \text (phage/mL)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{r_{max}}{\varphi _0} - 10^{10}$$\end{document} a.Bacteria lysis \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{0} \text (h)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.25 - 3$$\end{document} [9, 22, 24, 25] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{min} \text (h)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 - \frac{\tau _0}{2}$$\end{document} N/APhage release \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{mean} \text (phage/bact)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10 - 500$$\end{document} [24–28] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{\beta } \text (/h)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - 10$$\end{document} [14] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D \text (h)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 - \tau _0$$\end{document} b. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta \tau _{50} \text (h)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D - 3$$\end{document} c.Decays \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _{U} \text (/h)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.1 - 0.25$$\end{document} [29] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _{I} \text (/h)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.1 - 0.25$$\end{document} d. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _{P} \text (/h)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.05 - 0.15$$\end{document} [30, 31]Initial inoculums \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U_{init} \text (bact/mL)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^3 - K_C \cdot (1 - \frac{\delta _U}{r_{max}})$$\end{document} e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_{init} \text (phage/mL)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^3 - 10^9$$\end{document} N/Aa. Chosen to have a stronger maximal bactericidal effect than natural bacterial growth. b. Assuming phages have time to start replicating before lysing the bacteria. c. Taken in the range of allowed latent periods. d. Assuming similar range as uninfected bacteria decay. e. Smaller than maximal bacterial proliferation

Numerical implementation

The models were implemented in Python 3.1, using standard Python libraries (math, numpy, scipy, matplotlib, pandas). The random seed was set to 120659 for replicability purposes.

The differential equations were solved numerically with the forward Euler method adapted to distributed delay differential equations. The Runge-Kutta 4 method was also implemented but showed no noticeable convergence benefit compared to the Euler method for the full DDDE model (see appendix). The integrals were computed as convolution sums using the history of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi \cdot U \cdot P$$\end{document} (see appendix). Convergence was considered achieved if dividing the step size by two induced less than a 1% difference for both the bacterial load and the viral load profiles. When convergence failed during the Monte Carlo process, a new set of parameters was generated.

To correct for truncation of the normal distribution while preserving the desired mean and standard deviation, the Newton-Raphson method was used to find appropriate parameters.

Results

Distributions of \documentclass[12pt]{minimal}

			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$\beta$$\end{document} and \documentclass[12pt]{minimal}
			\usepackage{amsmath}
			\usepackage{wasysym} 
			\usepackage{amsfonts} 
			\usepackage{amssymb} 
			\usepackage{amsbsy}
			\usepackage{mathrsfs}
			\usepackage{upgreek}
			\setlength{\oddsidemargin}{-69pt}
			\begin{document}$$\varphi$$\end{document} do not impact population dynamics

The analytical study demonstrated (see appendix) that the distributions of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi$$\end{document} can be simplified in the differential equations, and only their respective means \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _0(\tau )$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi (t)$$\end{document} impact population dynamics. The differential Eq. 7 could thus be simplified to:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{dU}{dt}(t)&= r(t) \cdot U(t) \nonumber \\&\qquad \quad - \varphi _0(t) \cdot U(t) \cdot P(t) - \delta _{U} \cdot U(t)\nonumber \\ \frac{dI}{dt}(t)&= \varphi _0(t) \cdot U(t) \cdot P(t) \\&\qquad \quad - \int _{\tau = 0}^{t} e^{- \delta _I \tau } \cdot \varphi _0(t-\tau ) \cdot U(t-\tau ) \cdot P(t-\tau ) \cdot f_{\tau }(\tau ) d\tau - \delta _I \cdot I(t) \nonumber \\ \frac{dP}{dt}(t)&= \int _{\tau = 0}^{t} \beta _{0}(\tau ) \cdot e^{- \delta _I \tau } \cdot \varphi _0(t-\tau ) \cdot U(t-\tau ) \cdot P(t-\tau ) \cdot f_{\tau }(\tau ) d\tau \nonumber \\&\qquad \qquad \qquad \qquad \qquad \qquad \quad - \varphi _0(t) \cdot (U(t) + I(t)) \cdot P(t)- \delta _P \cdot P(t) \nonumber \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _{0}(t)$$\end{document} is the mean adsorption rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi$$\end{document} at time t, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{0}(\tau )$$\end{document} is the mean burst size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} for a given latent period \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} . \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _0(t)$$\end{document} is assumed to follow Eq. 4 [13], and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _0(\tau )$$\end{document} to follow Eq. 5 [14]. Hereafter, Eq. 10 is used as the canonical expression for the full DDDE model.

The shape of the distribution has negligible impact on population dynamics

Modest differences were observed visually in the bacterial and viral dynamics simulated with the lognormal, normal and gamma distributions (Fig. A3), even for wide distributions ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$rsd_{\tau } = 50\%$$\end{document} , higher end of observed latent period variability [17])).

This observation was quantitatively verified over 500 simulation scenarios (Fig. 5). For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$rsd_{\tau } = 50\%$$\end{document} , it represented less than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$6.3\%$$\end{document} median difference on bacterial load and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8.7\%$$\end{document} median difference on viral load. Differences were observed to be larger between the log-normal and normal distributions, which is consistent with those two distribution profiles being the furthest apart (Fig. A3). In all cases, the shape of the distribution had a generally limited impact on population dynamics.Fig. 5. Differences in bacterial and viral load predictions between log-normally, normally and gamma distributed latent period for equal distribution width, compiled over 500 random sets of parameters

Approximations with fixed latent period

Since the shape of the distribution was found to affect reasonably little the population dynamics, a distribution shape can be chosen arbitrarily to estimate the error induced from fixing the latent period. Here, the log-normal distribution was chosen to report the results.

Poor concordance was observed between the reference with distributed latent period and fixing the latent period to early observed lysis (Fig. 6 and A4). Even at moderate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$rsd_{\tau } = 10\%$$\end{document} , it already caused median bacterial and viral load prediction errors of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$18.4\%$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$18.4\%$$\end{document} respectively, and these errors increased rapidly with greater variability (Fig. 6B). Therefore, fixing the latent period \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} to the usually reported time to first lysis introduced a high bias in the model predictions (consistent with findings from [9, 11]).Fig. 6. Error on bacterial and viral load predictions induced by fixing the latent period to early lysis time, to the median latent period, or to the mean latent period, compiled over 500 random sets of parameters. The reference used is a log-normally distributed latent period

Fixing the latent period to the median or the mean of the distribution aligned much more closely with the reference with log-normally distributed latent period. The median error on both bacterial and viral load predictions remained limited below \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$12.0\%$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$15.3\%$$\end{document} , respectively, for a latent period variability up to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$rsd_{\tau } = 20\%$$\end{document} , which is a typically reported range of observed variability [4, 16]. This approximation did not hold at higher latent period variability, with up to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$29.8\%$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$37.8\%$$\end{document} median error respectively for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$rsd_{\tau } = 50\%$$\end{document} when fixing the latent period to its mean value.

TC models do not perform better than DDE models

As for the approximation with the TC model, it predicted smooth dynamics similar to the full DDDE model (Fig. A5). For low latent period variability up to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$rsd_{\tau } = 20\%$$\end{document} , this approximation resulted in reasonably low error with median of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$13.5\%$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$22.1\%$$\end{document} for bacterial and viral load predictions respectively (Fig. 7). This approximation however also broke at high variability with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$53.1\%$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$63.1\%$$\end{document} respective median error at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$rsd_{\tau } = 50\%$$\end{document} . Overall, the TC model did not perform better than the Fixed Mean DDE model, which demonstrated the highest performance among the DDE models.Fig. 7. Error on bacterial and viral load predictions of the DDE and TC models relative to the full DDDE model, compiled over 500 random sets of parameters. The reference used is a log-normally distributed latent period

Discussion

We successfully proposed a novel model structure that can account for cell-level variability as well as intracellular phage replication dynamics (Fig. 3 & Eq. 10). The random variabilities on the burst size and the adsorption rate were shown to not impact population-level dynamics, thus only the variability on the latent period and its influence on the burst size remained in the final model proposed. We then used this model structure to quantify systematically said impact to help inform more precise model building for phage dynamics modelling. A wide latent period distribution (such as the one observed for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document} X174 [17]) was shown to greatly influence the population-level dynamics. As for lower variability, more traditional approaches, with fixed latent period or with transit compartments, were able to predict bacterial and viral load profiles close to the ones obtained with the full model.

Since the developed model structure is a distributed delay differential equation (DDDE) model, it is computationally intensive (7s to 9s of runtime, see appendix). For narrower latent period distributions, both the fixed delay differential equation (DDE) and the transit compartments (TC) models can be close approximations of the full model. The TC model is a common usual ODE structure and can be implemented in any ODE software, while offering better run times than the full DDDE model (0.1s to 1s). However, it requires one more parameter than the DDE model, the number of transit compartments is sometimes hard to estimate, and for very narrow latent period ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$rsd_\tau \le 10\%$$\end{document} ), this number of transit compartments can be very high ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N = round(\frac{1}{rsd_\tau ^2}) \ge 100$$\end{document} [15]) resulting in potentially slower run times ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ge$$\end{document} 1s). As for the DDE model, it overall performed better than the TC model, provided that the latent period should be fixed to its mean or median value rather than the time of first observed lysis, to prevent a systematic bias, while offering very efficient run times ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx$$\end{document} 0.1s). However, DDE equations are not supported in every ODE softwares and might require some additional work to be implemented. At higher variability, the TC model introduced higher prediction errors than the DDE model in the general case. However, the former can perform exceptionally well even at high latent period variability if the burst size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} varies little with the latent period \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _0(\tau ) \approx constant$$\end{document} , for instance if we can neglect the dynamics relative to the intracellular viral replication (see appendix). As such, a priori knowledge of expected latent period variability and intracellular dynamics can greatly help inform model selection. In practice, estimating phage parameters from population dynamics in the DDE and TC models may result in lower errors than previously described. However, this may introduce bias on the other parameters (burst size, adsorption rate) to compensate for the latent period misspecification. An extension of this work could thus explore how well those two models can fit profiles simulated with the full DDDE model, and whether such models remain reliable for extrapolation beyond the range used for the parameter estimation.

The approach proposed here remains theoretical, based on a semi-mechanistic structural model. Despite using model parameters in typically observed ranges, the parameter combinations simulated might not always correspond to real bacteriophages, and we only compared a few usual distributions for the distributed parameters. Given available experimental data [4, 16], distributions were not expected to differ greatly from a normal distribution. As such only normal, log-normal and gamma distributions were considered in this work. Yet, other distributions, especially less-standard distributions, such as bimodal distributions, might result in different trends. The results discussed here should thus be further validated through experimental studies. While precise measurement of cell-level variability can require labor-intensive single-cell experiments and cannot be easily applied to wider pharmacological applications, population-level studies can still provide valuable insights into the underlying distribution, as demonstrated through single-step growth curves. Tracking bacterial load and viral load over a single round of infection could allow us to identify latent period distribution and the correlation between burst size and latent period, and assess the concordance of the developed model with real profiles.

This work can also provide insight on model selection in different contexts. For simulation purposes, run time is often a less critical criterion and the full DDDE model may be preferred since it is expected to capture the dynamics more mechanistically. Measuring cell-level heterogeneity, to properly inform the full DDDE model, can however remain challenging experimentally, and different assumptions and translation work from other phages might be required. For estimation purposes, both the DDE and the TC model offer a much more practical run time. Given that phages generally exhibit low variability in latent period [4, 16], it may be advisable to first try using the DDE model, which requires the least parameters, in absence of prior knowledge regarding the variability. In that case, it is preferable to either directly estimate a latent period from population dynamics, if feasible, or to fix it based on an experimentally determined median or mean value. A TC model can also be tested, especially in software with advanced traditional ODE solvers that might not be compatible with DDE models. In this case, the number of transit compartments can either be fixed to an arbitrary value higher than N = 10, to account for the fact that typical rsd is less than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$30\%$$\end{document} , or estimated. Conversely, if experiments show evidence of strong latent period variability, i.e., a gradual single-step growth curve, or if the previous two models fail to capture the profile accurately, more accurate distributed approaches may still be necessary. Some approaches relying on partial differential equations can for instance be used to allow for a thorough description of aging dynamics of the infected bacteria [32, 33], but those are highly computationally intensive and difficult to implement in some modeling softwares often used in pharmacometrics, such as NONMEM. The full DDDE model with intracellular replication dynamics could offer a good in-between, that efficiently captures the infection age of infected cells to track cell-level dynamics. We are also currently working on an ODE approximation of this approach to further facilitate integration into any ODE solver, improve significantly run time, and allow for a potentially broader use in the field.

In conclusion, the classical fixed delay differential equation model should thus be used in first intention for phage modelling. However, if prior knowledge of latent period variability is available or if the fixed delay approach fails, a distributed approach is warranted and the full DDDE model may be the preferred option in most cases.

Supplementary Information

Below is the link to the electronic supplementary material.Supplementary file 1 (zip 18 KB)