# Phage-antibiotic synergy inhibited by temperate and chronic virus
competition

Kylie J. Landa  Department of Mathematics, Statistics, and Computer Science,
St. Olaf College, Northfield, Minnesota 55057, USA  Lauren M. Mossman
Department of Mathematics, Statistics, and Computer Science, St. Olaf College,
Northfield, Minnesota 55057, USA  Rachel J. Whitaker  Department of
Microbiology, University of Illinois at Urbana-Champaign, Urbana, Illinois
61801, USA  Carl R. Woese Institute for Genomic Biology, University of
Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA  Zoi Rapti
Department of Mathematics, University of Illinois at Urbana-Champaign, Urbana,
Illinois 61801, USA  Carl R. Woese Institute for Genomic Biology, University
of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA  Sara M. Clifton
<EMAIL_ADDRESS>Department of Mathematics, Statistics, and Computer
Science, St. Olaf College, Northfield, Minnesota 55057, USA

###### Abstract

As antibiotic resistance grows more frequent for common bacterial infections,
alternative treatment strategies such as phage therapy have become more widely
studied in the medical field. While many studies have explored the efficacy of
antibiotics, phage therapy, or synergistic combinations of phages and
antibiotics, the impact of virus competition on the efficacy of antibiotic
treatment has not yet been considered. Here, we model the synergy between
antibiotics and two viral types, temperate and chronic, in controlling
bacterial infections. We demonstrate that while combinations of antibiotic and
temperate viruses exhibit synergy, competition between temperate and chronic
viruses inhibits bacterial control with antibiotics. In fact, our model
reveals that antibiotic treatment may counterintuitively _increase_ the
bacterial load when a large fraction of the bacteria develop antibiotic-
resistance.

bacteria, bacteriophage, phage, latent, lytic, infection, recovery,
resistance, Pseudomonas aeruginosa, cystic fibrosis, population dynamics,
mathematical model

##  I Introduction

The rising number of antibiotic resistant bacteria is causing an increase in
the cost of treatment, higher mortality rates, and a longer period of sickness
duse2013 . For instance, infections by the multi-drug resistant pathogen
Pseudomonas aeruginosa, which is categorized as a serious threat by the CDC
CDC2019 , are frequently diagnosed among those with compromised immune systems
to such an extent that it is the most common cause of death among patients
with cystic fibrosis hancock2000antibiotic .

As antibiotic resistant bacteria are becoming more prevalent, new strategies
are needed to combat these bacterial infections. Recent experiments show that
antibiotics and phages may synergistically target bacterial populations when
used in tandem chaudhry2017synergy ; easwaran2020 . While some studies have
revealed that two phage types can synergistically fight bacterial infections
chaudhry2017synergy , we demonstrate that certain virus types that are
naturally present in human hosts may inhibit bacterial infection control with
antibiotics. Because the viruses compete for access to host cells and once the
cells are infected the viruses exclude one another, thus restricting the
sensitive population for the other virus.

The two viral lifestyles studied in this paper are temperate and chronic. When
chronic viruses infect a bacterium, the host cell produces and releases new
viruses without lysing itself weinbauer2004ecology . Chronic viruses also have
a latent (lysogenic) cycle in which viral genetic material is incorporated
into the bacterium’s genome. This genetic code is then transmitted to daughter
cells when the bacterium replicates lwoff1953lysogeny  . Unlike chronic
viruses, temperate viruses have a lytic and a latent cycle. During the latent
cycle, the virus remains dormant in the bacterium until induced to replicate.
During the lytic cycle, the bacterium produces viruses that burst out of the
cell, causing the host to die weinbauer2004ecology . This phenomenon of lysis
can be triggered by an SOS response induced in the bacterium cell
weinbauer2004ecology .

The SOS response is an inducible repair system found in bacteria such as P.
aeruginosa. Stress from a bacterium’s environment such as UV light, sublethal
doses of antibiotics, and radiation can precipitate this evolved response
fothergill2011effect ; rokney2008host . When certain classes of antibiotics,
such as quinolones, are introduced into the bacteria’s environment, they have
the potential to damage the bacterial DNA which can then induce the SOS
response power1992 . Once the SOS response has been activated, bacteria may
start to replicate previously latent viruses, potentially leading to cell
death Brives2020 .

This paper uses a mathematical model to investigate the effects of the two
competing virus types, temperate and chronic, on the effectiveness of various
antibiotic dosing frequencies. In most cases, we find that as the frequency of
antibiotic dosages increases, the bacterial population decreases even if the
bacteria are known to be antibiotic resistant. However, our model reveals that
antibiotic treatment counterintuitively _increases_ the bacterial load when a
large fraction of the bacteria develop antibiotic-resistance and temperate and
chronic phages are present in the system.

##  II Model

We begin with a model of temperate and chronic virus competition
clifton2019temperate , which assumes that two virus types, temperate viruses
($V_{T}$) with lytic and latent lytic stages and chronic viruses ($V_{C}$)
with productive and latent chronic stages, may infect a single strain of
bacteria that is initially susceptible ($S$) to both viral types. See Figure 1
for a model overview.

Figure 1: Compartmental model of bacteria-phage system with temperate and
chronic phages. Solid arrows represent bacteria transitions between infection
states, and dotted arrow represent virus production or loss. The green boxes
illustrate actions which may occur only when environmental stress, such as
sublethal antibiotics, are introduced to the system.

The total bacterial population is composed of susceptible, lytic, latent
temperate, productive, and latent chronic bacteria. We assume that the total
bacterial population grows logistically to a carrying capacity $K$
zwietering1990modeling . In the model system, we rescale all population
densities to be proportional to the carrying capacity, and we rescale the
timescale such that all rates are relative to the susceptible bacteria’s
growth rate.

The susceptible bacteria grow at an intrinsic rate $r_{S}$ and can become
infected by chronic or temperate viruses. The susceptible bacteria are
infected by the temperate viruses at a rate $\eta_{T}$. When infected, there
is a $f_{T}$ probability that the bacterium enters a latent lytic state in
which viruses are not being produced but the viral DNA is integrated into the
bacterium’s genome and passed on to later generations harper2014bacteriophages
. The latent lytic bacteria reproduce at a rate $r_{T}$ and recover from
infection at a rate $\gamma_{T}$.

Alternatively, there is a $(1-f_{T})$ probability that the bacteria will
instead enter the lytic lifestyle following infection. The bacteria hosting
phages in a lytic lifestyle die due to lysis (bursting of cell) at a rate
$\delta$. When the bacteria lyse, they produces $\beta_{T}$ temperate viruses
per burst bacterium. It is important to note that in both cases of recovery,
the bacteria become susceptible to phage infection once again. The bacteria
can recover from infection at a rate $\gamma_{I}$.

Similar to temperate infections, the susceptible bacteria are infected by the
chronic viruses at a rate $\eta_{C}$. When infected, there is a $f_{C}$
probability that the bacteria will enter the latent chronic state in which the
bacteria can reproduce at a rate $r_{C}$ and recover from infection at a rate
$\gamma_{C}$.

Alternatively, there is a $(1-f_{C})$ probability that the bacteria will
instead become productive upon infection and grow at a rate $r_{P}$ while
producing viruses at a rate $\beta_{C}$. In contrast to lytic infections,
productive infections do not cause cell death. Productive bacteria recover
from infection at a rate $\gamma_{P}$.

Outside the bacterial cells, free temperate viruses and free chronic viruses
are adsorbed by susceptible bacteria at a rate $\eta_{T}$ and $\eta_{C}$,
respectively, and degrade naturally at a rate $\mu_{T}$ and $\mu_{C}$,
respectively heldal1991production . Free temperate viruses are produced when
the lytic bacteria burst at a rate $\beta_{T}\delta$, and free chronic viruses
are generated by productive bacteria at a rate $\beta_{C}$.

Once a bacterium is infected, we assume it will exclude both superinfection by
the same viral type and cross infection by viruses of the other type
de2017pseudomonas . Although we assume infection rates for temperate and
chronic viruses are constant, the infection rates may depend on the bacterial
population; many relevant bacteria form biofilms at high density that protect
the population from infection harper2014bacteriophages . For simplicity, we
have also assumed that lysogen frequencies are constant, but some studies have
demonstrated that bacterial density may impact lysogeny rates
hargreaves2014does ; silpe2018host . These simplifications are necessary for
analytic tractability.

###  II.1 Antibiotics

We introduce antibiotics into the system at times $t_{i}$; specifically these
are the times at which the antibiotics become bioavailable, which depends on
the mode of delivery. We assume that system stress spikes at times $t_{i}$ and
decays exponentially, consistent with typical antibiotic metabolism in the
human system naber1973pharmacokinetics ; bax1989pharmacokinetics ;
fish1997clinical . The antibiotic concentration within the system is then:

$\displaystyle a(t)=A\,\sum_{i=1}^{N}H(t-t_{i})\,e^{-k(t-t_{i})},$ (1)

where $t$ is the current time, $t_{i}$ are the antibiotic dose times, $A$ is
the maximum concentration of one antibiotic dose, $N$ is the total number of
antibiotic doses, $H$ is the Heaviside function, and $k$ is the decay rate of
antibiotics in the system.

When the system is stressed by antibiotics, there are three major responses.
First, antibiotics can induce the production of latent temperate or chronic
viruses at a rate $\epsilon\,a(t)$ fothergill2011effect ; ptashne1986genetic ;
nanda2015impact . This response occurs with bacteria that are antibiotic
resistant or antibiotic susceptible fisher2017persistent ;
monack2004persistent ; stewart2001antibiotic ; james2012differential ;
redgrave2014fluoroquinolone ; valencia2017ciprofloxacin ;
brazas2005ciprofloxacin . However, not all classes of antibiotics induce virus
production, so this model is only applicable to antibiotics (such as
quinolones) known to do so zhang2000quinolone ; comeau2007phage . To simplify
the model, we ignore spontaneous induction because it is a rare phenomenon
nanda2015impact ; garro1974relationship ; Cortes546275 .

Second, antibiotic resistant bacteria are not killed directly by the
antibiotics fisher2017persistent ; monack2004persistent ;
stewart2001antibiotic , but antibiotic susceptible bacteria are killed at a
rate proportional to the antibiotic concentration: $\kappa\,a(t)$
levin2010population .

Finally, studies show that productive bacteria stressed by sublethal
antibiotics increase viral production while decreasing cell reproduction
hagens2006augmentation ; secor2015filamentous . We incorporate the increased
viral production with a Hollings-like functional response; when there is no
antibiotic stress in the system ($a=0$), viruses are produced at a rate
$\beta_{C}$, and as the concentration of antibiotics increases, the viral
production rate saturates at $\beta_{\mathrm{max}}$ (see Figure 2a). The
functional form of the stress-dependent production rate is

$\displaystyle\beta_{C}(a)=\beta_{C}+\frac{a}{h_{\beta}+a}(\beta_{\mathrm{max}}-\beta_{C}),$
(2)

where $a$ is the time-dependent antibiotic concentration, $h_{b}$ is the
concentration at which the production rate is halfway between the minimum and
maximum, and $\beta_{max}$ is the maximum production rate when stress is
maximum.

Figure 2: Functional forms of chronic virus production rates (left) and
productive bacteria growth rates (right). With no system stress (no
antibiotics), the chronic phage production rate is $\beta_{C}$ and the
bacterial growth rate is $r_{P}$. As stress increases, the phage production
rate saturates to $\beta_{\mathrm{max}}$ and the bacterial growth rate
approaches 0.

We also incorporate the decreased cell reproduction when the concentration of
antibiotics increases using a Hollings-like functional response. When there is
no stress to the system, the cellular reproduction rate is $r_{p}$. As the
amount of stress (antibiotics) increases, the reproduction rate approaches 0
(See Figure 2b). The functional form of the stress-dependent cellular
reproduction rate is

$\displaystyle r_{P}(a)=r_{P}\left(1-\frac{a}{h_{r}+a}\right),$ (3)

where $a$ is the time-dependent antibiotic concentration, and $h_{r}$ is the
concentration at which the reproduction rate is halfway between the minimum
and maximum.

The dynamical system capturing the described model behavior is

$\displaystyle\dot{S}=\underbrace{S\left(1-N\right)}_{\text{growth}}-\underbrace{\eta_{T}SV_{T}-\eta_{C}SV_{C}}_{\text{infection}}+\underbrace{\gamma_{T}L_{T}+\gamma_{I}I_{T}+\gamma_{P}P_{C}+\gamma_{C}L_{C}}_{\text{recovery}}-\underbrace{\kappa
a(t)S}_{\begin{subarray}{c}\text{antibiotic}\\\ \text{death}\end{subarray}}$
(4)
$\displaystyle\dot{I}_{T}=\underbrace{\eta_{T}(1-f_{T})SV_{T}}_{\text{infection}}-\underbrace{\delta
I_{T}}_{\text{lysis}}-\underbrace{\gamma_{I}I_{T}}_{\text{recovery}}+\underbrace{\epsilon
a(t)L_{T}}_{\text{induction}}-\underbrace{\kappa
a(t)I_{T}}_{\begin{subarray}{c}\text{antibiotic}\\\
\text{death}\end{subarray}}$ (5)
$\displaystyle\dot{L}_{T}=\underbrace{r_{T}L_{T}\left(1-N\right)}_{\text{growth}}+\underbrace{\eta_{T}f_{T}SV_{T}}_{\text{infection}}-\underbrace{\gamma_{T}L_{T}}_{\text{recovery}}-\underbrace{\epsilon
a(t)L_{T}}_{\text{induction}}-\underbrace{\kappa
a(t)L_{T}}_{\begin{subarray}{c}\text{antibiotic}\\\
\text{death}\end{subarray}}$ (6)
$\displaystyle\dot{P}_{C}=\underbrace{r_{P}(a)P_{C}\left(1-N\right)}_{\text{growth}}+\underbrace{(1-f_{C})\eta_{C}SV_{C}}_{\text{infection}}-\underbrace{\gamma_{P}P_{C}}_{\text{recovery}}+\underbrace{\epsilon
a(t)L_{C}}_{\text{induction}}-\underbrace{\kappa
a(t)P_{C}}_{\begin{subarray}{c}\text{antibiotic}\\\
\text{death}\end{subarray}}$ (7)
$\displaystyle\dot{L}_{C}=\underbrace{r_{C}L_{C}\left(1-N\right)}_{\text{growth}}+\underbrace{f_{C}\eta_{C}SV_{C}}_{\text{infection}}-\underbrace{\gamma_{C}L_{C}}_{\text{recovery}}-\underbrace{\epsilon
a(t)L_{C}}_{\text{induction}}-\underbrace{\kappa
a(t)L_{C}}_{\begin{subarray}{c}\text{antibiotic}\\\
\text{death}\end{subarray}}$ (8)
$\displaystyle\dot{V}_{T}=\underbrace{\beta_{T}\delta
I_{T}}_{\text{burst}}-\underbrace{\eta_{T}SV_{T}}_{\text{adsorption}}-\underbrace{\mu_{T}V_{T}}_{\text{degradation}}$
(9)
$\displaystyle\dot{V}_{C}=\underbrace{\beta_{C}(a)P_{C}}_{\text{production}}-\underbrace{\eta_{C}SV_{C}}_{\text{adsorption}}-\underbrace{\mu_{C}V_{C}}_{\text{degradation}}$
(10)

All model variables are listed in Table 1 and all model parameters are listed
in Table 2.

Table 1: Description of model variables in bacteria-virus system (4-10). In the system, all population densities are rescaled to be proportional to the carrying capacity and the timescale is rescaled such that all rates are relative to the susceptible bacteria’s growth rate. See supplemental materials for rescaling process. Variable |  Meaning |  Units
---|---|---
$S$ |  susceptible bacteria as a proportion of carrying capacity |  unitless
$I_{T}$ |  lytic bacteria as a proportion of carrying capacity |  unitless
$L_{T}$ |  latent lytic bacteria as a proportion of carrying capacity |  unitless
$P_{C}$ |  productive bacteria as a proportion of carrying capacity |  unitless
$L_{C}$ |  latent chronic bacteria as a proportion of carrying capacity |  unitless
$N$ |  all bacteria ($S+I_{T}+L_{T}+P_{C}+L_{C}$) |  unitless
$V_{T}$ |  ratio of free temperate viruses to bacteria carrying capacity |  $\displaystyle\frac{\text{PFU}}{\text{CFU}}$
$V_{C}$ |  ratio of free chronic viruses to bacteria carrying capacity |  $\displaystyle\frac{\text{PFU}}{\text{CFU}}$
$t$ |  time rescaled by intrinsic growth rate |  unitless
Table 2: Description of model parameters in bacteria-virus system (4-10). See
supplemental materials for parameter selection process and model rescaling
process.

Parameter |  Meaning |  Rangea |  Baseline |  Sources
---|---|---|---|---
$r_{T},r_{C}$ |  growth rates of (respectively) latent lytic and latent chronic bacteria, relative to the susceptible bacteria growth rateb |  [0.5, 3]c |  1 |  shapiro2016evolution
$r_{P}(a)$ |  growth rate of productive bacteria, a decreasing function of the antibiotic concentration $a$ |  [0.5, 3]c |  1 |  shapiro2016evolution
$\eta_{T},\eta_{C}$ |  infection rate of (respectively) temperate and chronic viruses |  [0.38, 14.7]e |  1 |  sinha2017silico
$\gamma_{T},\gamma_{P},\gamma_{C}$ |  recovery rates of (respectively) latent lytic, productive, and latent chronic bacteria |  [0,1000]f |  0.67g |  brussow2004phages ; horvath2010crispr
$\gamma_{I}$ |  recovery rates of lytic bacteria |  [0,1000] |  0h |  brussow2004phages ; horvath2010crispr
$\delta$ |  rate at which lytic infection leads to bursting |  [1.5, 7.8]i |  4 |  yu2017characterization ; el2015isolation
$f_{T}$ |  lysogen frequency for temperate viruses |  [0, 0.9] |  0.01 |  calendar2006bacteriophages ; oppenheim2007new ; volkova2014modeling
$f_{C}$ |  lysogen frequency for chronic viruses |  [0, 0.9] |  0g |  calendar2006bacteriophages ; oppenheim2007new ; volkova2014modeling
$\beta_{T}$ |  burst size for bacteria infected with $V_{T}$ |  [10, 1000] |  100 |  yu2017characterization ; el2015isolation ; latino2014novel ; schrader1997bacteriophage ; ceyssens2010molecular ; garbe2011sequencing ; you2002effects
$\beta_{C}(a)$ |  viral production rate for bacteria infected with $V_{C}$, an increasing function of antibiotic concentration $a$ |  [5, 200] |  20 |  clifton2019modeling
$\mu_{T},\mu_{C}$ |  degradation rate of (respectively) free temperate viruses and free chronic viruses |  [0.9, 3.6]j |  1 |  heldal1991production
$\kappa$ |  death rate per concentration of antibiotics |  [0, 3.5] |  1 |  spalding2018mathematical ; grillon2016comparative
$\epsilon$ |  lysis rate per concentration of antibiotics |  [0, 2] |  1 |  spalding2018mathematical ; grillon2016comparative
$A$ |  maximum concentration of one standard antibiotic dose ($\mu$g/mL) |  [0, 10]k |  1 |  grillon2016 ; fong1986
$k$ |  metabolic decay rate of antibiotic within the system |  [1e-3,0.6]l |  0.3 |  spalding2018mathematical ; zhanel2006review ; wingender1984pharmacokinetics
$h_{\beta},h_{r}$ |  stress level at which growth rate is half the maximum |  |  1 |

aall parameter ranges are taken for the human pathogens P. aeruginosa or E.
coli and their viruses, unless otherwise noted. bgrowth rate is approximately
5.1e-3 min-1 for P. aeruginosa grown in vitro, but is highly variable in human
hosts. cestimates based on E. coli and M13 phage. dstable bacterial density in
human hosts is highly variable; a study of viable P. aeruginosa densities in
sputum of 12 patients with cystic fibrosis not undergoing treatment ranged
from 5.3e3 CFU/mL to 1.8e11 CFU/mL. eestimates based on E. coli and $\lambda$
phage. fa wide range of recovery rates has been found; some proviruses are
viable over evolutionary timescales and some proviruses are inactivated nearly
instantly by CRISPR systems. gestimated from viral steady state density (see
Results section). hselected to be 0 to simplify model analysis; allowing
$\gamma_{I}=\gamma_{T}$ produces qualitatively similar results, so the
increased model complexity is not justified. ilow estimate is for PAXYB1 phage
and PAO1 host, high estimate is for $\varphi$PSZ1 phage and PAO1 host. jlow
estimate is for viruses extracted from Raunefjorden, high estimate is for
viruses extracted from Bergen Harbor (strains unknown). kblood serum
concentration of ciprofloxacin and levofloxacin. lantibiotic is levofloxacin.

##  III Results

Because patients infected with dangerous bacterial infections, particularly P.
aeruginosa, are typically treated with antibiotics at the time of bacterial
detection hoiby2017diagnosis ; hoiby2015escmid , we investigate the effects of
antibiotic administration on the bacteria-phage ecosystem. We assume that the
antibiotics trigger an SOS response in the bacterium which induces the
production of latent temperate and chronic viruses fothergill2011effect ;
ptashne1986genetic ; nanda2015impact . Not all antibiotics induce phage via
the SOS response fothergill2011effect , so we focus only on the types of
antibiotics known to do so (e.g., quinolones like levofloxacin and
ciprofloxacin) zhang2000quinolone ; comeau2007phage . We assume that even
antibiotic-resistant bacteria induce viruses in the presence of antibiotics,
which has been demonstrated for several classes of antibiotics
redgrave2014fluoroquinolone ; valencia2017ciprofloxacin ;
brazas2005ciprofloxacin ; james2012differential ; fothergill2011effect .

We also assume that as the concentration of antibiotics increases, the
production rate of the chronic viruses increases and the reproduction rate of
productive bacteria decreases hagens2006augmentation ; secor2015filamentous .
This assumption is critical to the main results presented here, so we
thoroughly explore the effect of chronic virus production as a function of
antibiotic administration frequency in the Sensitivity Analysis section.

Figure 3 shows the time-averaged bacterial population level as a function of
the frequency of antibiotic administration. Antibiotic susceptible bacterial
populations are represented with solid lines, and antibiotic resistant
bacterial populations are represented with dashed lines.

Figure 3: Simulation of the effect of antibiotic administration frequency on
the time-averaged bacterial population when infected with various
bacteriophages in the human host. Time is scaled by the growth rate of P.
aeruginosa ($5.1$e$-3$ min-1), so 0.136 in unitless frequency units is
approximately one dose per 24 hours. In the system, all population densities
are rescaled to be proportional to the carrying capacity. Thus, the “placebo”
condition (no antibiotic treatment) results in the bacterial population
reaching carrying capacity (1 in unitless density).

###  III.1 Antibiotic Susceptible

Figure 3 illustrates that no matter the antibiotic frequency, the bacterial
population is largest when not infected with viruses. When the SOS response is
triggered in bacteria without bacteriophages, no lysis from rapid
bacteriophage reproduction occurs. Instead, the bacteria can only die from
antibiotics. Alternatively, when bacteriophages are present, the bacterial
population decreases at a faster rate. When the frequency of antibiotics is
below 0.12, the presence of temperate and chronic viruses has a larger
negative impact on the bacterial population than when only chronic viruses are
present. However, once the antibiotic frequency is larger than 0.12, the
bacterial populations become equivalent when temperate and chronic and only
chronic are present. When no antibiotics are present, the chronic bacterial
population reaches carrying capacity because the bacteria are not lysing or
being killed by antibiotics. However, when chronic and temperate viruses are
present with no antibiotics, the temperate viruses can still lyse the
bacteria. This accounts for the difference in bacterial population when
frequency of antibiotic administration is zero. The two populations converge
when the antibiotic frequency is greater than 0.12 because chronic viruses
out-compete temperate viruses due to the stress-induced chronic virus
production rate increase.

Antibiotics have the largest effect on the bacterial population infected with
only temperate viruses. When infected with temperate viruses, the bacterial
population severely decreases as frequency of antibiotics increases because
latent bacteria are being induced by antibiotics into the lytic lifestyle. The
apparent noise in population levels when the frequency is greater than 0.08 is
likely due to chaos in the system. Studies show that chaos is commonly seen in
periodically forced predator-prey models tang2003 ; taylor2013 . Further
exploration of this phenomenon is beyond the scope of this paper.

###  III.2 Antibiotic Resistant

When no bacteriophages are present in the system or when only chronic viruses
are present in the system, the antibiotic-resistant bacterial populations
remain at carrying capacity as the frequency of antibiotics increases. The
bacterial population reaches carrying capacity when no viruses are present
because there exists no mechanism (antibiotic or virus) to kill the bacteria.
The bacterial population also reaches carrying capacity when only chronic
viruses are present because chronic viruses do not kill their host in order to
reproduce.

The behaviors of the antibiotic-resistant bacterial population change when
both viruses are in the system and when only temperate viruses are in the
system. The bacterial population infected by only temperate viruses decreases
as the frequency of antibiotics is increased. Decreasing at an increasing
rate, we see that the bacterial population levels off at 0.11. This is because
as the frequency of antibiotics increases, more latent temperate viruses are
induced into the lytic cycle. The dynamics of the antibiotic resistant
bacteria infected with temperate bacteriophages are similar to when the
bacteria are antibiotic susceptible, however when the bacteria are antibiotic
resistant, it requires a higher frequency of antibiotics to significantly
decrease the bacterial population.

Contrary to expectations, as antibiotic frequency increases, the antibiotic-
resistant bacterial population infected with temperate and chronic viruses
increases. At a frequency of roughly 0.13, the bacterial population infected
with temperate and chronic viruses reaches carrying capacity. While
antibiotics are not killing the antibiotic resistant bacteria through the
intended mechanisms, they cause the induction of the SOS response in bacteria.
This induction causes the production of chronic viruses and lysis of temperate
infected bacteria. At higher antibiotic frequencies, all of the temperate
infected bacteria lyse. The only bacteria that are left are the bacteria that
have been infected by chronic viruses. Chronic viruses, having out-competed
temperate viruses (a natural predator of bacteria), leave bacteria
invulnerable to control. This phenomenon has clinical implications in that,
contrary to the intended outcome, antibiotics can counterintuitively increase
the bacterial population.

###  III.3 Sensitivity Analysis

While we have used P. aeruginosa infections as a case study, this model is
theoretically applicable to any ecosystem with bacteria, temperate viruses,
and chronic viruses. Therefore, we conduct a sensitivity and uncertainty
analysis to understand the impact of each model parameter on the clinical
outcome of interest: bacteria population size. Because the relationships
between the parameters and the bacterial population size are generally
nonlinear but monotonic, we use Latin Hypercube Sampling (LHS) of parameter
space and partial rank correlation coefficients (PRCC). Following Marino et
al. marino2008methodology , the sensitivity analysis indicates that the
parameters $\kappa$, $T$, $\beta_{max}$, and $\eta$ have a significant
($p<0.05$) impact on the bacterial population (Figure 4). The antibiotic death
rate ($\kappa$), has the largest magnitude of PRCC, indicating that a small
increase in $\kappa$ will result in the largest decrease of the bacterial
population. The period of antibiotic administration ($T$), the infection rate
$\eta$, and the maximum chronic virus production rate $\beta_{max}$ also have
a negative correlation with the bacterial population size.

For any particular bacteria-phage-antibiotic system, the parameters
$\kappa,T$, and $\eta$ should be well-characterised. It is less likely that
$\beta_{max}$ will be well-characterized because chronic viruses have not been
studied as extensively as lytic and temperate viruses. We wish to understand
the impact of $\beta_{max}$ on the bacterial population when the population is
antibiotic resistant and infected with both temperate and chronic viruses;
this scenario is the most common for a patient infected with P. aeruginosa
james2015lytic ; burgener2019 ; geller2011levofloxacin . Figure 5 shows that
as $\beta_{max}$ decreases, the slope of bacterial population as a function
frequency of antibiotic decreases until it ultimately becomes negative around
$\beta_{max}=70$. This shows that virus competition is only detrimental to
infection control if chronic virus production is sufficiently increased by
antibiotic stress.

Figure 4: Sensitivity of parameters on average bacteria population. The
sensitivity analyses uses Latin hypercube sampling (LHS) of parameter space
and partial rank correlation coefficients (PRCC) marino2008methodology . All
parameter values are taken near the baselines in Table 2. Initially,
$S(0)=1e-3,V_{T}(0)=V_{C}(0)=1e-7$. 500 simulations were completed. Note that
asterisks indicate significance. See the supplemental material for technical
details. Figure 5: Impact of $\beta_{max}$ on resistant bacterial population
when temperate and chronic bacteriophages are present. As $\beta_{max}$
decreases, the effect of antibiotic frequency on bacterial population flips
from a positive slope to a negative slope.

##  IV Discussion

The implications of our model can theoretically be applied to any bacterial
infection with naturally present temperate and chronic viruses, but the
results presented in this paper focus on treatment of P. aeruginosa infections
in patients with cystic fibrosis. Clinical trials of antibiotics confirm that
both temperate and chronic viral lifestyles are present in the lungs of most
patients with cystic fibrosis james2015lytic ; burgener2019 , and around 60%
of P. aeruginosa in sputum samples are antibiotic resistant
geller2011levofloxacin . Despite the rise of antibiotic-resistant bacterial
infections, antibiotics remain a common treatment for patients with cystic
fibrosis proesmans2013 . Unfortunately, our model suggests that antibiotic
treatment may be counterproductive under the most likely treatment conditions.
These results highlight the need for personalized medicine approaches for
treating antibiotic-resistant bacterial infection which respect, or even
exploit, the underlying ecology within the patient.

###  IV.1 Limitations and Future Steps

Our results are limited by assumptions made in order to create a simple,
analytically tractable model. First, the model assumes that the antibiotics
present induce viral production through the SOS response. Not all antibiotics
cause an SOS response in bacteria. Thus, the model is applicable to
antibiotics known to trigger an SOS response (e.g., quinolone antibiotics like
levofloxacin and ciprofloxacin) fothergill2011effect . This is not a major
limitation as treatment of P. aeruginosa infections often include quinolone
antibiotics geller2011levofloxacin ; hodson1987oral .

The model also assumes mass action infection of bacteria by phage. Due to the
biological nature of bacteria, this assumption breaks down at large bacterial
densities. Bacteria like P. aeruginosa are known to form biofilms at high
densities, which effectively saturates the infection rate
harper2014bacteriophages . Previous studies have used a Hollings Type II
functional response to model the infection rate, and the results are not
qualitatively different clifton2019modeling . However, a more sophisticated
model would include biofilm formation in the infection assumptions.

The phenomenon of multiple infections, where more than one phage infects a
single bacterium at one time is excluded from our model. This is not a
significant limitation because when infected with bacteriophages, P.
aeruginosa produces superinfection exclusion proteins which prevent multiple
infections by the same phage type james2012differential ; heo2007genome .
Another limitation of our model is that we do not consider the severity of an
infection. The concept of multiplicity, where more viruses infecting the host
may increase the severity of the illness, is not considered because the
severity of infection is not included in our model.

From the sensitivity analysis, we know that parameters $\kappa$, $T$,
$\beta_{max}$, and $\eta$ have the largest impact on the accuracy of our
model. While $T$ is determined by the clinician and $\kappa$ and $\eta$ are
generally well-characterized for a broad range of phages and antibiotics,
$\beta_{max}$ has not been definitively determined experimentally. Figure 5
illustrates the impact of $\beta_{max}$ on antibiotic resistant bacteria
infected with temperate and chronic viruses. Figure 5 illustrates the
importance of having a precise measurement of $\beta_{max}$. Further
experiments are needed to determine a precise value of $\beta_{max}$, which
would allow the model to be more readily applicable in a clinical setting.

Finally, because our model does not incorporate evolutionary dynamics, such as
the development of antibiotic or phage resistance, our results are only
applicable on short time scales. Our model could serve as the base of a more
sophisticated model that includes both ecological and evolutionary dynamics.

##  V Conclusion

This model demonstrates the potential for an antagonistic effect of virus
competition on antibiotic treatment of resistant bacterial infections. Both
temperate and chronic virus lifestyles are naturally present in human hosts,
and antibiotic resistance is a growing concern. Because certain combinations
of phages and antibiotics can exhibit a synergistic effect when treating
bacterial infections, this work suggests that a more personalized medicine
approach may provide better clinical outcomes. Treatment plans tailored to the
underlying bacteria-phage ecology naturally present within human hosts are
critical to controlling bacterial infections, especially when the bacteria are
resistant to antibiotics.

##  VI Acknowledgements

This work was funded in part by the National Science Foundation grant
DMS-1815764 (ZR), the Cystic Fibrosis Foundation grant WHITAK16PO (RJW), and
an Allen Distinguished Investigator Award (RJW). The funders had no role in
study design, data collection and analysis, decision to publish, or
preparation of the manuscript.

The authors thank Jayadevi H. Chandrashekhar and George A. O’Toole for
discussions that informed biological aspects of this work. Thanks are also due
to Ted Kim and Karna Gowda for assistance with parameter selection. The
authors also thank Laura Suttenfield and Alan Collins for suggestions that
improved the clarity of the paper.

##  VII Data availability

All software (Matlab .m files) will be made publicly available via the
Illinois Data Bank.

##  VIII Competing interests

The authors declare no competing interests.

## References

  * (1) R. Laxminarayan, A. Duse, C. Wattal, A. K. Zaidi, H. Wertheim, N. Sumpradit, E. Vlieghe, G. L. Hara, I. M. Gould, H. Goossens, C. Greko, A. D. So, M. Bigdeli, G. Tomson, W. Woodhouse, E. Ombaka, A. Q. Peralta, F. Qamar, F. Mir, S. Kariuki, Z. A. Bhutta, A. Coates, R. Bergstrom, G. D. Wright, E. D. Brown, and C. Otto, “Antibiotic resistance - the need for global solutions,” The Lancet, vol. 13, no. 12, pp. 1057–1098, 2013.
  * (2) C. for Disease Control and Prevention, “2019 ar threats report,” 2019.
  * (3) R. E. Hancock and D. P. Speert, “Antibiotic resistance in pseudomonas aeruginosa: mechanisms and impact on treatment,” Drug resistance updates, vol. 3, no. 4, pp. 247–255, 2000.
  * (4) W. N. Chaudhry, J. Concepción-Acevedo, T. Park, S. Andleeb, J. J. Bull, and B. R. Levin, “Synergy and order effects of antibiotics and phages in killing pseudomonas aeruginosa biofilms,” PloS one, vol. 12, no. 1, p. e0168615, 2017.
  * (5) M. Easwaran, M. De Zoysa, and H.-J. Shin, “Application of phage therapy: Synergistic effect of phage ecsw and antibiotic combination towards antibiotic-resistant escherichia coli,” Transboundary and Emerging Diseases, 2020.
  * (6) M. G. Weinbauer, “Ecology of prokaryotic viruses,” FEMS microbiology reviews, vol. 28, no. 2, pp. 127–181, 2004.
  * (7) A. Lwoff, “Lysogeny,” Bacteriological Reviews, vol. 17, no. 4, p. 269, 1953\.
  * (8) J. L. Fothergill, E. Mowat, M. J. Walshaw, M. J. Ledson, C. E. James, and C. Winstanley, “Effect of antibiotic treatment on bacteriophage production by a cystic fibrosis epidemic strain of pseudomonas aeruginosa,” Antimicrobial agents and chemotherapy, vol. 55, no. 1, pp. 426–428, 2011.
  * (9) A. Rokney, O. Kobiler, A. Amir, D. L. Court, J. Stavans, S. Adhya, and A. B. Oppenheim, “Host responses influence on the induction of lambda prophage,” Molecular microbiology, vol. 68, no. 1, pp. 29–36, 2008.
  * (10) E. Power and I. Phillips, “Induction of the sos gene (umuc) by 4-quinolone antibacterial drugs,” Journal of Medical Microbiology, vol. 36, no. 2, 1992\.
  * (11) C. Brives and J. Pourraz, “Phage therapy as a potential solution in the fight against amr: obstacles and possible futures,” Palgrave Commun, vol. 6, no. 100, 2020.
  * (12) S. M. Clifton, R. J. Whitaker, and Z. Rapti, “Temperate and chronic virus competition leads to low lysogen frequency,” arXiv preprint arXiv:1911.07233, 2019.
  * (13) M. Zwietering, I. Jongenburger, F. Rombouts, and K. Van’t Riet, “Modeling of the bacterial growth curve,” Applied and environmental microbiology, vol. 56, no. 6, pp. 1875–1881, 1990.
  * (14) D. R. Harper, H. M. Parracho, J. Walker, R. Sharp, G. Hughes, M. Werthén, S. Lehman, and S. Morales, “Bacteriophages and biofilms,” Antibiotics, vol. 3, no. 3, pp. 270–284, 2014.
  * (15) M. Heldal and G. Bratbak, “Production and decay of viruses in aquatic environments,” Marine Ecology Progress Series, pp. 205–212, 1991.
  * (16) J. De Smet, H. Hendrix, B. G. Blasdel, K. Danis-Wlodarczyk, and R. Lavigne, “Pseudomonas predators: Understanding and exploiting phage–host interactions,” Nature Reviews Microbiology, vol. 15, no. 9, p. 517, 2017\.
  * (17) K. R. Hargreaves, A. M. Kropinski, and M. R. Clokie, “What does the talking?: quorum sensing signalling genes discovered in a bacteriophage genome,” PloS one, vol. 9, no. 1, p. e85131, 2014.
  * (18) J. E. Silpe and B. L. Bassler, “A host-produced quorum-sensing autoinducer controls a phage lysis-lysogeny decision,” Cell, 2018.
  * (19) K. G. Naber, S. R. Westenfelder, and P. O. Madsen, “Pharmacokinetics of the aminoglycoside antibiotic tobramycin in humans,” Antimicrobial agents and chemotherapy, vol. 3, no. 4, pp. 469–473, 1973.
  * (20) R. Bax, W. Bastain, A. Featherstone, D. Wilkinson, and M. Hutchison, “The pharmacokinetics of meropenem in volunteers,” Journal of antimicrobial chemotherapy, vol. 24, no. suppl_A, pp. 311–320, 1989.
  * (21) D. N. Fish and A. T. Chow, “The clinical pharmacokinetics of levofloxacin,” Clinical pharmacokinetics, vol. 32, no. 2, pp. 101–119, 1997.
  * (22) M. Ptashne, A genetic switch: Gene control and phage. lambda.  Palo Alto, CA (US); Blackwell Scientific Publications, 1986.
  * (23) A. M. Nanda, K. Thormann, and J. Frunzke, “Impact of spontaneous prophage induction on the fitness of bacterial populations and host-microbe interactions,” Journal of bacteriology, vol. 197, no. 3, pp. 410–419, 2015\.
  * (24) R. A. Fisher, B. Gollan, and S. Helaine, “Persistent bacterial infections and persister cells,” Nature Reviews Microbiology, vol. 15, no. 8, p. 453, 2017\.
  * (25) D. M. Monack, A. Mueller, and S. Falkow, “Persistent bacterial infections: the interface of the pathogen and the host immune system,” Nature Reviews Microbiology, vol. 2, no. 9, p. 747, 2004.
  * (26) P. S. Stewart and J. W. Costerton, “Antibiotic resistance of bacteria in biofilms,” The lancet, vol. 358, no. 9276, pp. 135–138, 2001.
  * (27) C. E. James, J. L. Fothergill, H. Kalwij, A. J. Hall, J. Cottell, M. A. Brockhurst, and C. Winstanley, “Differential infection properties of three inducible prophages from an epidemic strain of pseudomonas aeruginosa,” BMC microbiology, vol. 12, no. 1, p. 216, 2012.
  * (28) L. S. Redgrave, S. B. Sutton, M. A. Webber, and L. J. Piddock, “Fluoroquinolone resistance: mechanisms, impact on bacteria, and role in evolutionary success,” Trends in microbiology, vol. 22, no. 8, pp. 438–445, 2014.
  * (29) E. Y. Valencia, F. Esposito, B. Spira, J. Blázquez, and R. S. Galhardo, “Ciprofloxacin-mediated mutagenesis is suppressed by subinhibitory concentrations of amikacin in pseudomonas aeruginosa,” Antimicrobial agents and chemotherapy, vol. 61, no. 3, pp. e02107–16, 2017.
  * (30) M. D. Brazas and R. E. Hancock, “Ciprofloxacin induction of a susceptibility determinant in pseudomonas aeruginosa,” Antimicrobial agents and chemotherapy, vol. 49, no. 8, pp. 3222–3227, 2005.
  * (31) X. Zhang, A. D. McDaniel, L. E. Wolf, G. T. Keusch, M. K. Waldor, and D. W. Acheson, “Quinolone antibiotics induce shiga toxin-encoding bacteriophages, toxin production, and death in mice,” The Journal of infectious diseases, vol. 181, no. 2, pp. 664–670, 2000.
  * (32) A. M. Comeau, F. Tétart, S. N. Trojet, M.-F. Prere, and H. Krisch, “Phage-antibiotic synergy (pas): $\beta$-lactam and quinolone antibiotics stimulate virulent phage growth,” PLoS One, vol. 2, no. 8, p. e799, 2007\.
  * (33) A. J. Garro and M.-F. Law, “Relationship between lysogeny, spontaneous induction, and transformation efficiencies in bacillus subtilis,” Journal of bacteriology, vol. 120, no. 3, pp. 1256–1259, 1974.
  * (34) M. G. Cortes, J. Krog, and G. Balazsi, “Optimality of the spontaneous prophage induction rate,” bioRxiv, 2019.
  * (35) B. R. Levin and K. I. Udekwu, “Population dynamics of antibiotic treatment: a mathematical model and hypotheses for time-kill and continuous-culture experiments,” Antimicrobial agents and chemotherapy, vol. 54, no. 8, pp. 3414–3426, 2010.
  * (36) S. Hagens, A. Habel, and U. Bläsi, “Augmentation of the antimicrobial efficacy of antibiotics by filamentous phage,” Microbial Drug Resistance, vol. 12, no. 3, pp. 164–168, 2006.
  * (37) P. R. Secor, J. M. Sweere, L. A. Michaels, A. V. Malkovskiy, D. Lazzareschi, E. Katznelson, J. Rajadas, M. E. Birnbaum, A. Arrigoni, K. R. Braun, et al., “Filamentous bacteriophage promote biofilm assembly and function,” Cell host & microbe, vol. 18, no. 5, pp. 549–559, 2015.
  * (38) J. W. Shapiro, E. S. Williams, and P. E. Turner, “Evolution of parasitism and mutualism between filamentous phage m13 and escherichia coli,” PeerJ, vol. 4, p. e2060, 2016.
  * (39) V. Sinha, A. Goyal, S. L. Svenningsen, S. Semsey, and S. Krishna, “In silico evolution of lysis-lysogeny strategies reproduces observed lysogeny propensities in temperate bacteriophages,” Frontiers in microbiology, vol. 8, p. 1386, 2017.
  * (40) H. Brüssow, C. Canchaya, and W.-D. Hardt, “Phages and the evolution of bacterial pathogens: from genomic rearrangements to lysogenic conversion,” Microbiol. Mol. Biol. Rev., vol. 68, no. 3, pp. 560–602, 2004.
  * (41) P. Horvath and R. Barrangou, “Crispr/cas, the immune system of bacteria and archaea,” Science, vol. 327, no. 5962, pp. 167–170, 2010.
  * (42) X. Yu, Y. Xu, Y. Gu, Y. Zhu, and X. Liu, “Characterization and genomic study of Òphikmv-likeÓ phage paxyb1 infecting pseudomonas aeruginosa,” Scientific reports, vol. 7, no. 1, p. 13068, 2017.
  * (43) G. El Didamony, A. Askora, and A. A. Shehata, “Isolation and characterization of t7-like lytic bacteriophages infecting multidrug resistant pseudomonas aeruginosa isolated from egypt,” Current microbiology, vol. 70, no. 6, pp. 786–791, 2015.
  * (44) R. Calendar, The bacteriophages.  Oxford University Press on Demand, 2006.
  * (45) A. B. Oppenheim, S. L. Adhya, et al., “A new look at bacteriophage $\lambda$ genetic networks,” Journal of bacteriology, vol. 189, no. 2, pp. 298–304, 2007.
  * (46) V. V. Volkova, Z. Lu, T. Besser, and Y. T. Gröhn, “Modeling the infection dynamics of bacteriophages in enteric escherichia coli: estimating the contribution of transduction to antimicrobial gene spread,” Appl. Environ. Microbiol., vol. 80, no. 14, pp. 4350–4362, 2014.
  * (47) L. Latino, C. Essoh, Y. Blouin, H. V. Thien, and C. Pourcel, “A novel pseudomonas aeruginosa bacteriophage, ab31, a chimera formed from temperate phage paju2 and p. putida lytic phage af: characteristics and mechanism of bacterial resistance,” PloS one, vol. 9, no. 4, p. e93777, 2014.
  * (48) H. S. Schrader, J. O. Schrader, J. J. Walker, T. A. Wolf, K. W. Nickerson, and T. A. Kokjohn, “Bacteriophage infection and multiplication occur in pseudomonas aeruginosa starved for 5 years,” Canadian journal of microbiology, vol. 43, no. 12, pp. 1157–1163, 1997.
  * (49) P.-J. Ceyssens, A. Brabban, L. Rogge, M. S. Lewis, D. Pickard, D. Goulding, G. Dougan, J.-P. Noben, A. Kropinski, E. Kutter, et al., “Molecular and physiological analysis of three pseudomonas aeruginosa phages belonging to the Òn4-like virusesÓ,” Virology, vol. 405, no. 1, pp. 26–30, 2010\.
  * (50) J. Garbe, B. Bunk, M. Rohde, and M. Schobert, “Sequencing and characterization of pseudomonas aeruginosa phage jg004,” BMC microbiology, vol. 11, no. 1, p. 102, 2011.
  * (51) L. You, P. F. Suthers, and J. Yin, “Effects of escherichia coli physiology on growth of phage t7 in vivo and in silico,” Journal of bacteriology, vol. 184, no. 7, pp. 1888–1894, 2002.
  * (52) S. M. Clifton, T. Kim, J. H. Chandrashekhar, G. A. O’Toole, Z. Rapti, and R. J. Whitaker, “Lying in wait: Modeling the control of bacterial infections via antibiotic-induced proviruses,” MSystems, vol. 4, no. 5, pp. e00221–19, 2019.
  * (53) C. Spalding, E. Keen, D. J. Smith, A.-M. Krachler, and S. Jabbari, “Mathematical modelling of the antibiotic-induced morphological transition of pseudomonas aeruginosa,” PLoS computational biology, vol. 14, no. 2, p. e1006012, 2018.
  * (54) A. Grillon, F. Schramm, M. Kleinberg, and F. Jehl, “Comparative activity of ciprofloxacin, levofloxacin and moxifloxacin against klebsiella pneumoniae, pseudomonas aeruginosa and stenotrophomonas maltophilia assessed by minimum inhibitory concentrations and time-kill studies,” PloS one, vol. 11, no. 6, p. e0156690, 2016.
  * (55) A. Grillon, F. Schramm, A. Vandenbroucke, M. Simbul, and V. Rahm, “Comparative activity of ciprofloxacin, levofloxacin and moxifloxacin against klebsiella pneumoniae, pseudomonas aeruginosa and stenotrophomonas maltophilia assessed by minimum inhibitory concentrations and time-kill studies,” PloS one, vol. 11, no. 6, 2016.
  * (56) I. Fong, W. Ledbetter, M. Kleinberg, and F. Jehl, “Ciprofloxacin concentrations in bone and muscle after oral dosing,” Plos One, vol. 29, no. 3, p. 405–408, 1986.
  * (57) G. G. Zhanel, S. Fontaine, H. Adam, K. Schurek, M. Mayer, A. M. Noreddin, A. S. Gin, E. Rubinstein, and D. J. Hoban, “A review of new fluoroquinolones,” Treatments in respiratory medicine, vol. 5, no. 6, pp. 437–465, 2006.
  * (58) W. Wingender, K.-H. Graefe, W. Gau, D. Förster, D. Beermann, and P. Schacht, “Pharmacokinetics of ciprofloxacin after oral and intravenous administration in healthy volunteers,” European Journal of Clinical Microbiology, vol. 3, no. 4, pp. 355–359, 1984.
  * (59) N. Høiby, T. Bjarnsholt, C. Moser, P. Ø. Jensen, M. Kolpen, T. Qvist, K. Aanæs, T. Pressler, M. Skov, and O. Ciofu, “Diagnosis of biofilm infections in cystic fibrosis patients,” Apmis, vol. 125, no. 4, pp. 339–343, 2017.
  * (60) N. Høiby, T. Bjarnsholt, C. Moser, G. Bassi, T. Coenye, G. Donelli, L. Hall-Stoodley, V. Hola, C. Imbert, K. Kirketerp-Møller, et al., “Escmid? guideline for the diagnosis and treatment of biofilm infections 2014,” Clinical microbiology and infection, vol. 21, pp. S1–S25, 2015\.
  * (61) S. Tang and L. Chen, “Quasiperiodic solutions and chaos in a periodically forced predator-prey model with age structure for predator,” International Journal of Bifurcation and Chaos, vol. 13, no. 4, pp. 973–980, 2003.
  * (62) R. Taylor, J. Sherratt, and A. White, “Seasonal forcing and multi-year cycles in interacting populations: lessons from a predator-prey model,” Journal of mathematical biology, vol. 67, no. 6, pp. 1741–1764, 2013.
  * (63) S. Marino, I. B. Hogue, C. J. Ray, and D. E. Kirschner, “A methodology for performing global uncertainty and sensitivity analysis in systems biology,” Journal of theoretical biology, vol. 254, no. 1, pp. 178–196, 2008.
  * (64) C. E. James, E. V. Davies, J. L. Fothergill, M. J. Walshaw, C. M. Beale, M. A. Brockhurst, and C. Winstanley, “Lytic activity by temperate phages of pseudomonas aeruginosa in long-term cystic fibrosis chronic lung infections,” The ISME journal, vol. 9, no. 6, p. 1391, 2015.
  * (65) E. Burgener, J. Sweere, M. Bach, P. Secor, N. Haddock, L. Jennings, R. Marvig, H. K. Johansen, E. Rossi, X. Cao, L. Tian, L. Nedelec, S. Molin, P. Bollyky, and C. Milla, “Filamentous bacteriophages are associated with chronic pseudomonas lung infections and antibiotic resistance in cystic fibrosis,” Science Translational Medicine, vol. 11, no. 488, 2019.
  * (66) D. E. Geller, P. A. Flume, D. Staab, R. Fischer, J. S. Loutit, and D. J. Conrad, “Levofloxacin inhalation solution (mp-376) in patients with cystic fibrosis with pseudomonas aeruginosa,” American journal of respiratory and critical care medicine, vol. 183, no. 11, pp. 1510–1516, 2011.
  * (67) M. Proesmans, F. Vermeulen, L. Boulanger, J. Verhaegen, and K. De Boeck, “Comparison of two treatment regimens for eradication of pseudomonas aeruginosa infection in children with cystic fibrosis,” Journal of Cystic Fibrosis, vol. 12, no. 1, pp. 29–34, 2012.
  * (68) M. Hodson, R. Butland, C. Roberts, M. Smith, and J. Batten, “Oral ciprofloxacin compared with conventional intravenous treatment for pseudomonas aeruginosa infection in adults with cystic fibrosis,” The Lancet, vol. 329, no. 8527, pp. 235–237, 1987.
  * (69) Y.-J. Heo, I.-Y. Chung, K. B. Choi, G. W. Lau, and Y.-H. Cho, “Genome sequence comparison and superinfection between two related pseudomonas aeruginosa phages, d3112 and mp22,” Microbiology, vol. 153, no. 9, pp. 2885–2895, 2007.
  * (70) S. H. Kopf, A. L. Sessions, E. S. Cowley, C. Reyes, L. Van Sambeek, Y. Hu, V. J. Orphan, R. Kato, and D. K. Newman, “Trace incorporation of heavy water reveals slow and heterogeneous pathogen growth rates in cystic fibrosis sputum,” Proceedings of the National Academy of Sciences, vol. 113, no. 2, pp. E110–E116, 2016.
  * (71) F. A. Stressmann, G. B. Rogers, P. Marsh, A. K. Lilley, T. W. Daniels, M. P. Carroll, L. R. Hoffman, G. Jones, C. E. Allen, N. Patel, N. Forbes, B. Forbes, A. Tuck, and K. D. Bruce, “Does bacterial density in cystic fibrosis sputum increase prior to pulmonary exacerbation?,” Journal of Cystic Fibrosis, vol. 10, no. 5, pp. 357–365, 2011.
  * (72) G. S. Stent et al., “Molecular biology of bacterial viruses.,” Molecular biology of bacterial viruses., 1963.
  * (73) J. J. Barr, R. Auro, N. Sam-Soon, S. Kassegne, G. Peters, N. Bonilla, M. Hatay, S. Mourtada, B. Bailey, M. Youle, et al., “Subdiffusive motion of bacteriophage in mucosal surfaces increases the frequency of bacterial encounters,” Proceedings of the National Academy of Sciences, vol. 112, no. 44, pp. 13675–13680, 2015.
  * (74) T. F. Thingstad, S. Våge, J. E. Storesund, R.-A. Sandaa, and J. Giske, “A theoretical analysis of how strain-specific viruses can control microbial species diversity,” Proceedings of the National Academy of Sciences, vol. 111, no. 21, pp. 7813–7818, 2014.
  * (75) D. Kirschner, “Uncertainty and sensitivity functions and implementation,” 2008\.
  * (76) M. McKay, R. Beckman, and W. Conover, “Comparison of three methods for selecting values of input variables in the analysis of output from a computer cod,” Technometrics, vol. 21, no. 2, p. 239–245, 1979.
  * (77) A. Saltelli, “Making best use of model evaluations to compute sensitivity indices,” Computer physics communications, vol. 145, no. 2, pp. 280–297, 2002.

##  IX Appendix

###  IX.1 Model rescaling

The dimensional model system (4-10) below can be rescaled so that time and
bacterial density are unitless quantities. Although the viral density can also
be nondimensionalized to be unitless, we have elected (for the sake of
interpretability) to rescale viral density to have units of PFU/CFU.

$\displaystyle\dot{\tilde{S}}=\underbrace{\tilde{r}_{S}\tilde{S}\left(1-\frac{\tilde{N}}{\tilde{K}}\right)}_{\text{growth}}-\underbrace{\tilde{\eta}_{T}\tilde{S}\tilde{V}_{T}-\tilde{\eta}_{C}\tilde{S}\tilde{V}_{C}}_{\text{infection}}+\underbrace{\tilde{\gamma}_{T}\tilde{L}_{T}+\tilde{\gamma}_{I}\tilde{I}_{T}+\tilde{\gamma}_{P}\tilde{P}_{C}+\tilde{\gamma}_{C}\tilde{L}_{C}}_{\text{recovery}}-\underbrace{\tilde{\kappa}\tilde{a}(\tilde{t})\tilde{S}}_{\begin{subarray}{c}\text{antibiotic}\\\
\text{death}\end{subarray}}$ (11)
$\displaystyle\dot{\tilde{I}}_{T}=\underbrace{\tilde{\eta}_{T}(1-\tilde{f}_{T})\tilde{S}\tilde{V}_{T}}_{\text{infection}}-\underbrace{\tilde{\delta}\tilde{I}_{T}}_{\text{lysis}}-\underbrace{\tilde{\gamma}_{I}\tilde{I}_{T}}_{\text{recovery}}+\underbrace{\tilde{\epsilon}\tilde{a}(\tilde{t})\tilde{L_{T}}}_{\text{induction}}-\underbrace{\tilde{\kappa}\tilde{a}(\tilde{t})\tilde{I_{T}}}_{\begin{subarray}{c}\text{antibiotic}\\\
\text{death}\end{subarray}}$ (12)
$\displaystyle\dot{\tilde{L}}_{T}=\underbrace{\tilde{r}_{T}\tilde{L}_{T}\left(1-\frac{\tilde{N}}{\tilde{K}}\right)}_{\text{growth}}+\underbrace{\tilde{\eta}_{T}\tilde{f}_{T}\tilde{S}\tilde{V}_{T}}_{\text{infection}}-\underbrace{\tilde{\gamma}_{T}\tilde{L}_{T}}_{\text{recovery}}-\underbrace{\tilde{\epsilon}\tilde{a}(\tilde{t})\tilde{L_{T}}}_{\text{induction}}-\underbrace{\tilde{\kappa}\tilde{a}(\tilde{t})\tilde{L_{T}}}_{\begin{subarray}{c}\text{antibiotic}\\\
\text{death}\end{subarray}}$ (13)
$\displaystyle\dot{\tilde{P}}_{C}=\underbrace{\tilde{r}_{P}\tilde{P}_{C}\left(1-\frac{\tilde{N}}{\tilde{K}}\right)}_{\text{growth}}+\underbrace{(1-\tilde{f}_{C})\tilde{\eta}_{C}\tilde{S}\tilde{V}_{C}}_{\text{infection}}-\underbrace{\tilde{\gamma}_{P}\tilde{P}_{C}}_{\text{recovery}}+\underbrace{\tilde{\epsilon}\tilde{a}(\tilde{t})\tilde{L_{C}}}_{\text{induction}}-\underbrace{\tilde{\kappa}\tilde{a}(\tilde{t})\tilde{P_{C}}}_{\begin{subarray}{c}\text{antibiotic}\\\
\text{death}\end{subarray}}$ (14)
$\displaystyle\dot{\tilde{L}}_{C}=\underbrace{\tilde{r}_{C}\tilde{L}_{C}\left(1-\frac{\tilde{N}}{\tilde{K}}\right)}_{\text{growth}}+\underbrace{\tilde{f}_{C}\tilde{\eta}_{C}\tilde{S}\tilde{V}_{C}}_{\text{infection}}-\underbrace{\tilde{\gamma}_{C}\tilde{L}_{C}}_{\text{recovery}}-\underbrace{\tilde{\epsilon}\tilde{a}(\tilde{t})\tilde{L_{C}}}_{\text{induction}}-\underbrace{\tilde{\kappa}\tilde{a}(\tilde{t})\tilde{L_{C}}}_{\begin{subarray}{c}\text{antibiotic}\\\
\text{death}\end{subarray}}$ (15)
$\displaystyle\dot{\tilde{V}}_{T}=\underbrace{\tilde{\beta}_{T}\tilde{\delta}\tilde{I}_{T}}_{\text{burst}}-\underbrace{\tilde{\eta}_{T}\tilde{S}\tilde{V}_{T}}_{\text{adsorption}}-\underbrace{\tilde{\mu}_{T}\tilde{V}_{T}}_{\text{degradation}}$
(16)
$\displaystyle\dot{\tilde{V}}_{C}=\underbrace{\tilde{\beta}_{C}\tilde{P}_{C}}_{\text{production}}-\underbrace{\tilde{\eta}_{C}\tilde{S}\tilde{V}_{C}}_{\text{adsorption}}-\underbrace{\tilde{\mu}_{C}\tilde{V}_{C}}_{\text{degradation}}$
(17)

Suppose equations (11-17) have bacterial density units of CFU/mL, viral
density units of PFU/mL, and time units of minutes; then equations (11-15)
have units of CFU/mL/min, and equations (16-17) have units of PFU/mL/min. We
will make the following substitutions into the system:

$\displaystyle\tilde{S}$ $\displaystyle=\tilde{K}S$
$\displaystyle\tilde{I}_{T}$ $\displaystyle=\tilde{K}I_{T}$
$\displaystyle\tilde{L}_{T}$ $\displaystyle=\tilde{K}L_{T}$
$\displaystyle\tilde{P}_{C}$ $\displaystyle=\tilde{K}P_{C}$
$\displaystyle\tilde{L}_{C}$ $\displaystyle=\tilde{K}L_{C}$
$\displaystyle\tilde{N}$ $\displaystyle=\tilde{K}N$
$\displaystyle\tilde{V}_{T}$ $\displaystyle=\tilde{K}V_{T}$
$\displaystyle\tilde{V}_{C}$ $\displaystyle=\tilde{K}V_{C}$
$\displaystyle\tilde{t}$ $\displaystyle=\frac{t}{\tilde{r}_{s}}$

where $\displaystyle S,\dots,N$ are unitless bacterial densities,
$\displaystyle V_{T},V_{C}$ have units of PFU/CFU (a virus to bacteria ratio),
and $t$ is a unitless time. An alternate substitution of
$\displaystyle\tilde{V}_{T}=\tilde{\beta}_{T}\tilde{K}V_{T}$ and
$\displaystyle\tilde{V}_{C}=\tilde{\beta}_{T}\tilde{K}V_{C}$ would have led to
a nondimensionalization of the system, but the nondimensionalized viral
densities would have been more challenging to interpret; therefore we eschew
this option.

We will illustrate the rescaling process for bacterial dynamics with equation
(11), and equations (12-15) are similar:

$\frac{\mathrm{d}\tilde{S}}{\mathrm{d}\tilde{t}}=\underbrace{\tilde{r}_{S}\tilde{S}\left(1-\frac{\tilde{N}}{\tilde{K}}\right)}_{\text{growth}}-\underbrace{\tilde{\eta}_{T}\tilde{S}\tilde{V}_{T}-\tilde{\eta}_{C}\tilde{S}\tilde{V}_{C}}_{\text{infection}}+\underbrace{\tilde{\gamma}_{T}\tilde{L}_{T}+\tilde{\gamma}_{I}\tilde{I}_{T}+\tilde{\gamma}_{P}\tilde{P}_{C}+\tilde{\gamma}_{C}\tilde{L}_{C}}_{\text{recovery}}-\underbrace{\tilde{\kappa}\tilde{a}(\tilde{t})\tilde{S}}_{\text{antibiotic
death}\\\ }$

After making the substitutions into equation (11), we get

$\displaystyle\frac{\mathrm{d}(\tilde{K}{S})}{\mathrm{d}({t}/\tilde{r}_{S})}=$
$\displaystyle\underbrace{\tilde{r}_{S}\tilde{K}{S}\left(1-\frac{\tilde{K}{N}}{\tilde{K}}\right)}_{\text{growth}}-\underbrace{\tilde{\eta}_{T}\tilde{K}{S}\tilde{K}{V}_{T}-\tilde{\eta}_{C}\tilde{K}{S}\tilde{K}{V}_{C}}_{\text{infection}}$
$\displaystyle+\underbrace{\tilde{\gamma}_{T}\tilde{K}{L}_{T}+\tilde{\gamma}_{I}\tilde{K}{I}_{T}+\tilde{\gamma}_{P}\tilde{K}{P}_{C}+\tilde{\gamma}_{C}\tilde{K}{L}_{C}}_{\text{recovery}}-\underbrace{\tilde{\kappa}\tilde{a}(\tilde{t})\tilde{K}S}_{\text{antibiotic
death}}$

We now divide the equation by $\tilde{K}\tilde{r}_{S}$:

$\frac{\mathrm{d}{S}}{\mathrm{d}{t}}=\underbrace{{S}\left(1-{N}\right)}_{\text{growth}}-\underbrace{\frac{\tilde{\eta}_{T}\tilde{K}}{\tilde{r}_{S}}{S}{V}_{T}-\frac{\tilde{\eta}_{C}\tilde{K}}{\tilde{r}_{S}}{S}{V}_{C}}_{\text{infection}}+\underbrace{\frac{\tilde{\gamma}_{T}}{\tilde{r}_{S}}{L}_{T}+\frac{\tilde{\gamma}_{I}}{\tilde{r}_{S}}{I}_{T}+\frac{\tilde{\gamma}_{P}}{\tilde{r}_{S}}{P}_{C}+\frac{\tilde{\gamma}_{C}}{\tilde{r}_{S}}{L}_{C}}_{\text{recovery}}-{\underbrace{\frac{\tilde{\kappa}\tilde{a}(\tilde{t})}{\tilde{r_{s}}}S}_{\text{antibiotic
death}}}$

The transformed parameters for growth, infection, and recovery are now
evident:

$\displaystyle{\eta}_{T}$
$\displaystyle=\frac{\tilde{\eta}_{T}\tilde{K}}{\tilde{r}_{S}}$
$\displaystyle{\eta}_{C}$
$\displaystyle=\frac{\tilde{\eta}_{C}\tilde{K}}{\tilde{r}_{S}}$
$\displaystyle{\gamma}_{T}$
$\displaystyle=\frac{\tilde{\gamma}_{T}}{\tilde{r}_{S}}$
$\displaystyle{\gamma}_{I}$
$\displaystyle=\frac{\tilde{\gamma}_{I}}{\tilde{r}_{S}}$
$\displaystyle{\gamma}_{P}$
$\displaystyle=\frac{\tilde{\gamma}_{P}}{\tilde{r}_{S}}$
$\displaystyle{\gamma}_{C}$
$\displaystyle=\frac{\tilde{\gamma}_{C}}{\tilde{r}_{S}}$

The transformed parameters for antibiotic death require more manipulation.
Because

$\tilde{a}(\tilde{t})=\tilde{A}\,\sum_{i=1}^{N}H(\tilde{t}-\tilde{t}_{i})\,e^{-\tilde{k}(\tilde{t}-\tilde{t}_{i})},$

we have the unitless parameter grouping111The argument of the Heaviside
function can be multiplied by any positive number without changing the output.

$\displaystyle\frac{\tilde{\kappa}\tilde{a}(\tilde{t})}{\tilde{r_{s}}}$
$\displaystyle=\frac{\tilde{\kappa}}{\tilde{r_{s}}}\tilde{A}\,\sum_{i=1}^{N}H(\tilde{r}_{S}\tilde{t}-\tilde{r}_{S}\tilde{t}_{i})\,e^{-\frac{\tilde{k}}{\tilde{r}_{S}}(\tilde{r}_{S}\tilde{t}-\tilde{r}_{S}\tilde{t}_{i})}$
$\displaystyle=\kappa\,A\,\sum_{i=1}^{N}H(t-t_{i})\,e^{-k(t-t_{i})}$
$\displaystyle=\kappa\,a(t),$

where the transformed parameters are

$\displaystyle\kappa$ $\displaystyle=\frac{\tilde{\kappa}}{\tilde{r_{s}}}$
$\displaystyle A$ $\displaystyle=\tilde{A}$ $\displaystyle k$
$\displaystyle=\frac{\tilde{k}}{\tilde{r_{s}}}$

With the previous substitutions, equation (11) reduces to

$\frac{\mathrm{d}{S}}{\mathrm{d}{t}}=\underbrace{{S}\left(1-{N}\right)}_{\text{growth}}-\underbrace{{\eta}_{T}{S}V_{T}-{\eta}_{C}{S}V_{C}}_{\text{infection}}+\underbrace{{\gamma}_{T}{L}_{T}+{\gamma}_{I}{I}_{T}+{\gamma}_{P}{P}_{C}+{\gamma}_{C}{L}_{C}}_{\text{recovery}}-\underbrace{\kappa
a(t)S}_{\begin{subarray}{c}\text{antibiotic}\\\ \text{death}\end{subarray}}$

Equations (12-15) are rescaled in a similar manner. Substituting the
transformations into equations (16-17) yields

$\displaystyle\frac{\mathrm{d}(\tilde{K}{V}_{T})}{\mathrm{d}({t}/\tilde{r}_{S})}$
$\displaystyle=\underbrace{\tilde{\beta}_{T}\tilde{\delta}\tilde{K}{I}_{T}}_{\text{burst}}-\underbrace{\tilde{\eta}_{T}\tilde{K}{S}\tilde{K}{V}_{T}}_{\text{adsorption}}-\underbrace{\tilde{\mu}_{T}\tilde{K}{V}_{T}}_{\text{degradation}}$
$\displaystyle\frac{\mathrm{d}(\tilde{K}{V}_{C})}{\mathrm{d}({t}/\tilde{r}_{S})}$
$\displaystyle=\underbrace{\tilde{\beta}_{C}\tilde{K}{P}_{C}}_{\text{production}}-\underbrace{\tilde{\eta}_{C}\tilde{K}{S}\tilde{K}{V}_{C}}_{\text{adsorption}}-\underbrace{\tilde{\mu}_{C}\tilde{K}{V}_{C}}_{\text{degradation}}$

Again dividing by $\tilde{K}\tilde{r}_{S}$, we get:

$\displaystyle\frac{\mathrm{d}{V}_{T}}{\mathrm{d}{t}}$
$\displaystyle=\underbrace{\tilde{\beta}_{T}\frac{\tilde{\delta}}{\tilde{r}_{S}}{I}_{T}}_{\text{burst}}-\underbrace{\frac{\tilde{\eta}_{T}\tilde{K}}{\tilde{r}_{S}}{S}{V}_{T}}_{\text{adsorption}}-\underbrace{\frac{\tilde{\mu}_{T}}{\tilde{r}_{S}}{V}_{T}}_{\text{degradation}}$
$\displaystyle\frac{\mathrm{d}{V}_{C}}{\mathrm{d}{t}}$
$\displaystyle=\underbrace{\frac{\tilde{\beta}_{C}}{\tilde{r}_{S}}{P}_{C}}_{\text{production}}-\underbrace{\frac{\tilde{\eta}_{C}\tilde{K}}{\tilde{r}_{S}}{S}{V}_{C}}_{\text{adsorption}}-\underbrace{\frac{\tilde{\mu}_{C}}{\tilde{r}_{S}}{V}_{C}}_{\text{degradation}}$

Adding to our list of transformed parameters

$\displaystyle{\beta}_{T}$ $\displaystyle=\tilde{\beta}_{T}$
$\displaystyle{\beta}_{C}$
$\displaystyle=\frac{\tilde{\beta}_{C}}{\tilde{r}_{S}}$
$\displaystyle{\mu}_{T}$ $\displaystyle=\frac{\tilde{\mu}_{T}}{\tilde{r}_{S}}$
$\displaystyle{\mu}_{C}$ $\displaystyle=\frac{\tilde{\mu}_{C}}{\tilde{r}_{S}}$
$\displaystyle{\delta}$ $\displaystyle=\frac{\tilde{\delta}}{\tilde{r}_{S}}$

the viral density equations reduce to

$\displaystyle\frac{\mathrm{d}{V}_{T}}{\mathrm{d}{t}}$
$\displaystyle=\underbrace{{\beta}_{T}{\delta}{I}_{T}}_{\text{burst}}-\underbrace{{\eta}_{T}{S}{V}_{T}}_{\text{adsorption}}-\underbrace{{\mu}_{T}{V}_{T}}_{\text{degradation}}$
$\displaystyle\frac{\mathrm{d}{V}_{C}}{\mathrm{d}{t}}$
$\displaystyle=\underbrace{{\beta}_{C}{P}_{C}}_{\text{production}}-\underbrace{{\eta}_{C}{S}{V}_{C}}_{\text{adsorption}}-\underbrace{{\mu}_{C}{V}_{C}}_{\text{degradation}}$

The full set of transformed parameters is listed in Table 3.

Table 3: Summary of model parameter groupings after rescaling of system (11-17). Parameter |  Units |  |  Associated parameter grouping |  Units
---|---|---|---|---
$\tilde{r}_{S}$ |  min-1 |  |  $\displaystyle{r}_{S}=\frac{\tilde{r}_{S}}{\tilde{r}_{S}}=1$ |  unitless
$\tilde{r}_{T},\tilde{r}_{P},\tilde{r}_{C}$ |  min-1 |  |  $\displaystyle{r}_{T}=\frac{\tilde{r}_{T}}{\tilde{r}_{S}},\,\,\,{r}_{P}=\frac{\tilde{r}_{P}}{\tilde{r}_{S}},\,\,\,{r}_{C}=\frac{\tilde{r}_{C}}{\tilde{r}_{S}}$ |  unitless
$\tilde{K}$ |  $\displaystyle\frac{\text{CFU}}{\text{mL}}$ |  |  $\displaystyle{K}=\frac{\tilde{K}}{\tilde{K}}=1$ |  unitless
$\tilde{\eta}_{T},\tilde{\eta}_{C}$ |  $\displaystyle\frac{\text{mL}}{\text{PFU min}}$ |  |  $\displaystyle{\eta}_{T}=\frac{\tilde{\eta}_{T}\tilde{K}}{\tilde{r}_{S}},\,\,\,{\eta}_{C}=\frac{\tilde{\eta}_{C}\tilde{K}}{\tilde{r}_{S}}$ |  $\displaystyle\frac{\text{CFU}}{\text{PFU}}$
$\tilde{\gamma}_{T},\tilde{\gamma}_{I},\tilde{\gamma}_{P},\tilde{\gamma}_{C}$ |  min-1 |  |  $\displaystyle{\gamma}_{T}=\frac{\tilde{\gamma}_{T}}{\tilde{r}_{S}},\,\,\,{\gamma}_{I}=\frac{\tilde{\gamma}_{I}}{\tilde{r}_{S}},\,\,\,{\gamma}_{P}=\frac{\tilde{\gamma}_{P}}{\tilde{r}_{S}},\,\,\,{\gamma}_{C}=\frac{\tilde{\gamma}_{C}}{\tilde{r}_{S}}$ |  unitless
$\tilde{\delta}$ |  min-1 |  |  $\displaystyle{\delta}=\frac{\tilde{\delta}}{\tilde{r}_{S}}$ |  unitless
$\tilde{f}_{T}$ |  unitless |  |  ${f}_{T}=\tilde{f}_{T}$ |  unitless
$\tilde{f}_{C}$ |  unitless |  |  ${f}_{C}=\tilde{f}_{C}$ |  unitless
$\tilde{A}$ |  $\mu$g/mL |  |  $A=\tilde{A}$ |  $\mu$g/mL
$\tilde{\kappa}$ |  $\displaystyle\frac{\text{mL}}{\mu\text{g min}}$ |  |  $\displaystyle\kappa=\frac{\tilde{\kappa}}{\tilde{r}_{S}}$ |  mL/$\mu$g
$\tilde{\epsilon}$ |  $\displaystyle\frac{\text{mL}}{\mu\text{g min}}$ |  |  $\displaystyle\epsilon=\frac{\tilde{\epsilon}}{\tilde{r}_{S}}$ |  mL/$\mu$g
$\tilde{k}$ |  min-1 |  |  $k=\frac{\tilde{k}}{\tilde{r}_{S}}$ |  unitless
$\tilde{\beta}_{T}$ |  $\displaystyle\frac{\text{PFU}}{\text{CFU}}$ |  |  ${\beta}_{T}=\tilde{\beta}_{T}$ |  $\displaystyle\frac{\text{PFU}}{\text{CFU}}$
$\tilde{\beta}_{C}$ |  $\displaystyle\frac{\text{PFU}}{\text{CFU min}}$ |  |  $\displaystyle{\beta}_{C}=\frac{\tilde{\beta}_{C}}{\tilde{r}_{S}}$ |  $\displaystyle\frac{\text{PFU}}{\text{CFU}}$
$\tilde{\mu}_{T},\tilde{\mu}_{C}$ |  min-1 |  |  $\displaystyle{\mu}_{T}=\frac{\tilde{\mu}_{T}}{\tilde{r}_{S}},\,\,\,{\mu}_{C}=\frac{\tilde{\mu}_{C}}{\tilde{r}_{S}}$ |  unitless

###  IX.2 Parameter selection

The growth rate $\tilde{r}_{S}$ for P. aeruginosa in vitro is approximately
5.1e-3 min-1 spalding2018mathematical , although P. aeruginosa growth is
highly variable in humans kopf2016trace . Therefore all rate parameters
provided in min-1 are divided by this rate in order to rescale (see Table 3).

The growth rates for latently infected bacteria, $\tilde{r}_{T}$ and
$\tilde{r}_{C}$, are assumed to be similar to uninfected bacteria. Our
estimate for the growth rate $\tilde{r}_{P}$ of productive bacteria is based
on E. coli infected with the filamentous phage M13 shapiro2016evolution ;
depending on the experimentally imposed fitness pressure, chronic viruses
appear capable of either increasing or decreasing host growth rates. In the
interest of simplicity, we assumed that chronic bacteria grow at a rate
similar to uninfected bacteria.

The carrying capacity $\tilde{K}$ of bacteria in a medium depends on the
environment. Even within the sputum of patients with cystic fibrosis, the
carrying capacity is difficult to estimate due to variability among patients.
One study of patients with cystic fibrosis found that the densities of viable
P. aeruginosa in sputum of 12 patients not undergoing treatment ranged from
5.3e3 CFU/mL to 1.8e11 CFU/mL stressmann2011does . Because of the wide range
of carrying capacities, we use the geometric mean of this range, 3e7 CFU/mL,
for illustrative purposes only. Because $\tilde{K}$ is challenging to
estimate, we use strategies that do not explicitly require $\tilde{K}$ to
estimate the infection rate.

The infection rate for E. coli and T4 phage in mucus (assuming mass action
infection) is known to be approximately 47e-10 mL/min per PFU
stent1963molecular ; barr2015subdiffusive . Infection in marine ecosystems
(also assuming mass action infection) is similar at around 24e-10 mL/min per
PFU stent1963molecular ; thingstad2014theoretical . Given our uncertainty in
the bacterial carrying capacity, we elect to use an infection rate within the
range given by Sinha et al. sinha2017silico ; the authors fit their mass
action infection model to time series population data that reached carrying
capacity. The authors present $\tilde{K}\tilde{\eta}\in[0.45,100]$ hr-1 and
$\tilde{r}_{S}\in[0.5,10]$ hr-1. Rescaling leads to a range of ${\eta}$
between $0.045$ and $200$, and our selected value of $\eta=1$ is near the
geometric mean of that range.

The phage production delay rate $\tilde{\delta}$ is estimated based on the
eclipse and rise phase of PAXYB1 and $\varphi$PSZ1 phage
yu2017characterization ; el2015isolation . The eclipse (latent) and rise phase
is 130 minutes total for PAXYB1 yu2017characterization  and 27 minutes total
for $\varphi$PSZ1 el2015isolation . The smallest (rescaled) delay rate is then
$1/130/$5.1e-3$=1.5$, and the largest is $1/27/$5.1e-3$=7.3$. We selected the
approximate average of this range, 4, to be the delay rate $\delta$.

The temperate lysogen frequency $f_{T}$ is estimated based on E. coli and
$\lambda$ phage oppenheim2007new : “It is known that a cell infected by one
phage predominantly follows the lytic default pathway (about 99% of the
time).”

The phage degradation rates $\tilde{\mu}_{T}$ and $\tilde{\mu}_{C}$ were
estimated based on the decay rates of phage in aquatic environments
heldal1991production . The decay rates ranged from 0.26 to 1.1 per hour. We
rescale by multiplying by 60 minutes per hour and the bacterial growth rate,
5.1e-3 per minute. The rescaled range of decay rates is then 0.9 to 3.6. We
selected a value of $\mu_{T}=\mu_{C}=1$ arbitrarily from this range.

The amplitude of stress $A$ is based on antibiotic concentrations ($\mu$g/mL)
in blood serum grillon2016 ; fong1986 . The antibiotics of interest in these
studies were ciprofloxacin and levofloxacin.

The metabolic decay rate of antibiotics $k$ is calculated based on the half-
life of levofloxacin zhanel2006review . The half-life of a standard dose of
levofloxacin within a human is approximately $T_{1/2}$ = 6.9 hours = 414
minutes; the metabolic decay rate is then $\tilde{k}=\ln{(2)}/T_{1/2}\approx
1.7$e-3 min-1; rescaling by the growth rate $\tilde{r}_{S}=5.1$e-3, we get
$k\approx 0.3$.

The values of $h_{\beta}$ and $h_{\gamma}$ were selected to be 1 for
simplicity. This means the antibiotic concentration for which the production
rate is halfway between the minimum and maximum is 1 $\mu$g/mL. Similarly, the
antibiotic concentration for which the growth rate is half the maximum is 1
$\mu$g/mL.

All other parameter baselines and ranges are educated guesses.

###  IX.3 Sensitivity analysis: technical details

In all analyses, we have used the initial conditions
$S(0)=1\mathrm{e}{-3},V_{T}(0)=V_{C}(0)=1\mathrm{e}{-7}$, following Sinha et
al. sinha2017silico .

The global uncertainty and sensitivity analysis was performed using the
methodology outlined in Marino et al. marino2008methodology . The base code
for the analysis is freely available at the author Denise Kirschner’s website
marinoCode .

In brief, the analysis uses Latin Hypercube Sampling (LHS) of parameter space
to simulate uncertainty in model parameters. LHS sampling requires fewer model
simulations than simple random sampling without introducing bias
mckay1979comparison . We used uniform sampling of each parameter about the
base values given in Table II. The ranges of the uniform samples are available
in our code.

We use Partial Rank Correlation Coefficients (PRCC) to test the sensitivity of
model outputs to parameter uncertainty because model outputs generally depend
monotonically on model inputs, but those relationships are not linear trends.
As noted by Marino et al. marino2008methodology , for linear trends we could
have used Pearson correlation coefficient (CC), partial correlation
coefficients (PCCs), or standardized regression coefficients (SRC). Had our
trends been non-monotonic, we would have used the Sobol method or one of its
many extensions saltelli2002making .

The displayed sensitivities in Fig 7 are the PRCCs for the model output $y$
(either minimum $T$ or minimum $\kappa$) and the model inputs $x_{j}$. As
described by Marino et al. marino2008methodology , partial rank correlation
characterizes the monotonic relationship between input $x_{j}$ and output $y$
after the effects on $y$ of the other inputs are removed. The values of PRCCs
fall between $-1$ and $1$, with $1$ indicating the strongest positive rank
correlation and $-1$ indicating the strongest negative rank correlation. The
significance indicates the probability that the rank correlation is zero
(i.e., large significance suggests that there is no relationship between
$x_{j}$ and $y$).


# PyFstat: a Python package for continuous gravitational-wave data analysis

David Keitel  Rodrigo Tenorio  Gregory Ashton  Reinhard Prix

††margin:  DOI:
https://doi.org/10.21105/joss.03000https://doi.org10.21105/joss.03000 Software
• https://github.com/openjournals/joss-
reviews/issues/3000https://doi.orgReview •
https://github.com/pyfstat/pyfstathttps://doi.orgRepository •
https://doi.org/10.5281/zenodo.4660591https://doi.orgArchive Editor:
https://danielskatz.org/https://doi.orgDaniel S. Katz
Reviewers: • https://github.com/RobertRoscahttps://doi.org@RobertRosca •
https://github.com/khanx169https://doi.org@khanx169 Submitted: 26 January 2021
Published: 06 April 2021 License
Authors of papers retain copyright and release the work under a Creative
Commons Attribution 4.0 International License
(http://creativecommons.org/licenses/by/4.0/https://doi.orgCC BY 4.0).

## Summary

Gravitational waves in the sensitivity band of ground-based detectors can be
emitted by a number of astrophysical sources, including not only binary
coalescences, but also individual spinning neutron stars. The most promising
signals from such sources, although as of 2020 not yet detected, are the long-
lasting, quasi-monochromatic ‘Continuous Waves’ (CWs). Many search methods
have been developed and applied on LIGO (Aasi et al. 2015) and Virgo (Acernese
et al. 2015) data. See Prix (2009), Riles (2017), and Sieniawska and Bejger
(2019) for reviews of the field.

The PyFstat package provides tools to perform a range of CW data analysis
tasks. It revolves around the $\mathcal{F}$-statistic, first introduced by
Jaranowski, Krolak, and Schutz (1998): a matched-filter detection statistic
for CW signals described by a set of frequency evolution parameters and
maximized over amplitude parameters. This has been one of the standard methods
for LIGO-Virgo CW searches for two decades. PyFstat is built on top of
established routines in LALSuite (LIGO Scientific Collaboration 2018) but
through its more modern Python interface it enables a flexible approach to
designing new search strategies.

Classes for various search strategies and target signals are contained in
three main submodules:

  * •

core: The basic wrappers to LALSuite’s $\mathcal{F}$-statistic algorithm. End-
users should rarely need to access these directly.

  * •

grid_based_searches: Classes to search over regular parameter-space grids.

  * •

mcmc_based_searches: Classes to cover promising parameter-space regions
through stochastic template placement with the Markov Chain Monte Carlo (MCMC)
sampler ptemcee (Vousden, Farr, and Mandel 2015).

Besides standard CWs from isolated neutron stars, PyFstat can also be used to
search for CWs from sources in binary systems (including the additional
orbital parameters), for CWs with a discontinuity at a pulsar glitch, and for
CW-like long-duration transient signals, e.g., from _after_ a pulsar glitch.
Specialized versions of both grid-based and MCMC-based search classes are
provided for these scenarios. Both fully-coherent and semi-coherent searches
(where the data is split into several segments for efficiency) are covered,
and an extension to the $\mathcal{F}$-statistic that is more robust against
single-detector noise artifacts (Keitel et al. 2014) is also supported. While
PyFstat’s grid-based searches do not compete with the sophisticated grid
setups and semi-coherent algorithms implemented in various LALSuite programs,
its main scientific use cases so far are for the MCMC exploration of
interesting parameter-space regions and for the long-duration transient case.

PyFstat was first introduced in Ashton and Prix (2018), which remains the main
reference for the MCMC-based analysis implemented in the package. The
extension to transient signals, which uses PyCUDA (Klöckner et al. 2012) for
speedup, is discussed in detail in Keitel and Ashton (2018), and the glitch-
robust search approaches in Ashton, Prix, and Jones (2018).

Additional helper classes, utility functions, and internals are included for
handling the common Short Fourier Transform (SFT) data format for LIGO data,
simulating artificial data with noise and signals in them, and plotting
results and diagnostics. Most of the underlying LALSuite functionality is
accessed through SWIG wrappings (Wette 2020) though for some parts, such as
the SFT handling, we still (as of the writing of this paper) call stand-alone
lalapps executables. Completing the backend migration to pure SWIG usage is
planned for the future.

The source of PyFstat is hosted on
https://github.com/PyFstat/PyFstat/https://doi.orgGitHub. The repository also
contains an automated test suite and a set of introductory example scripts.
Issues with the software can be submitted through GitHub and pull requests are
always welcome. PyFstat can be installed through pip, conda or docker
containers. Documentation in html and pdf formats is available from
https://readthedocs.org/projects/pyfstat/https://doi.orgreadthedocs.org and
installation instructions can be found there or in the
https://github.com/PyFstat/PyFstat/blob/master/README.mdhttps://doi.orgREADME
file. PyFstat is also listed in the Astrophysics Source Code Library as
https://ascl.net/2102.027https://doi.orgascl:2102.027.

## Statement of need

The sensitivity of searches for CWs and long-duration transient GWs is
generally limited by computational resources, as the required number of
matched-filter templates increases steeply for long observation times and wide
parameter spaces. The C-based LALSuite library (LIGO Scientific Collaboration
2018) contains many sophisticated search methods with a long development
history and high level of optimization, but is not very accessible for
researchers new to the field or for students; nor is it convenient for rapid
development and integration with modern technologies like GPUs or machine
learning. Hence, PyFstat serves a dual function of (i) making LALSuite CW
functionality more easily accessible through a Python interface, thus
facilitating the new user experience and, for developers, the exploratory
implementation of novel methods; and (ii) providing a set of production-ready
search classes for use cases not yet covered by LALSuite itself, most notably
for MCMC-based followup of promising candidates from wide-parameter-space
searches.

So far, PyFstat has been used for

  * •

the original proposal of MCMC followup for CW candidates (Ashton and Prix
2018);

  * •

developing glitch-robust CW search methods (Ashton, Prix, and Jones 2018);

  * •

speeding up long-transient searches with GPUs (Keitel and Ashton 2018);

  * •

followup of candidates from all-sky searches for CWs from sources in binary
systems, see Covas and Sintes (2020) and Abbott et al. (2021);

  * •

studying the impact of neutron star proper motions on CW searches (Covas
2020).

## Acknowledgements

We acknowledge contributions to the package from Karl Wette, Sylvia Zhu and
Dan Foreman-Mackey; as well as helpful suggestions by John T. Whelan, Luca
Rei, and the LIGO-Virgo-KAGRA Continuous Wave working group. D.K. and R.T. are
supported by European Union FEDER funds; the Spanish Ministerio de Ciencia,
Innovación y Universidades and Agencia Estatal de Investigación grants
PID2019-106416GB-I00/AEI/10.13039/501100011033, RED2018-102661-T,
RED2018-102573-E, FPA2017-90687-REDC, FPU 18/00694, and BEAGAL 18/00148
(cofinanced by the Universitat de les Illes Balears); the Comunitat Autonoma
de les Illes Balears through the Direcció General de Política Universitaria i
Recerca with funds from the Tourist Stay Tax Law ITS 2017-006 (PRD2018/24) and
the Conselleria de Fons Europeus, Universitat i Cultura; the Generalitat
Valenciana (PROMETEO/2019/071); and EU COST Actions CA18108, CA17137, CA16214,
and CA16104. This paper has been assigned document number LIGO-P2100008.

## References

Aasi, J., B. P. Abbott, R. Abbott, and others. 2015. “Advanced LIGO.” _Class.
Quant. Grav._ 32: 074001. https://doi.org/10.1088/0264-9381/32/7/074001.

Abbott, R., T. D. Abbott, S. Abraham, and others. 2021. “All-sky search in
early O3 LIGO data for continuous gravitational-wave signals from unknown
neutron stars in binary systems.” _Phys. Rev. D_ 103 (6): 064017\.
https://doi.org/10.1103/PhysRevD.103.064017.

Acernese, F., M. Agathos, K. Agatsuma, and others. 2015. “Advanced Virgo: a
second-generation interferometric gravitational wave detector.” _Class. Quant.
Grav._ 32 (2): 024001. https://doi.org/10.1088/0264-9381/32/2/024001.

Ashton, Gregory, and Reinhard Prix. 2018. “Hierarchical multistage MCMC
follow-up of continuous gravitational wave candidates.” _Phys. Rev. D_ 97
(10): 103020. https://doi.org/10.1103/PhysRevD.97.103020.

Ashton, Gregory, Reinhard Prix, and D.I. Jones. 2018. “A semicoherent glitch-
robust continuous-gravitational-wave search method.” _Phys. Rev. D_ 98 (6):
063011. https://doi.org/10.1103/PhysRevD.98.063011.

Covas, P. B. 2020. “Effects of proper motion of neutron stars on continuous
gravitational-wave searches.” _Mon. Not. Roy. Astron. Soc._ 500 (4): 5167–76.
https://doi.org/10.1093/mnras/staa3624.

Covas, P. B., and Alicia M. Sintes. 2020. “First all-sky search for continuous
gravitational-wave signals from unknown neutron stars in binary systems using
Advanced LIGO data.” _Phys. Rev. Lett._ 124 (19): 191102.
https://doi.org/10.1103/PhysRevLett.124.191102.

Jaranowski, Piotr, Andrzej Krolak, and Bernard F. Schutz. 1998. “Data analysis
of gravitational - wave signals from spinning neutron stars. 1. The Signal and
its detection.” _Phys. Rev. D_ 58: 063001.
https://doi.org/10.1103/PhysRevD.58.063001.

Keitel, David, and Gregory Ashton. 2018. “Faster search for long
gravitational-wave transients: GPU implementation of the transient
$\mathcal{F}$-statistic.” _Class. Quant. Grav._ 35 (20): 205003.
https://doi.org/10.1088/1361-6382/aade34.

Keitel, David, Reinhard Prix, Maria Alessandra Papa, Paola Leaci, and Maham
Siddiqi. 2014. “Search for continuous gravitational waves: Improving
robustness versus instrumental artifacts.” _Phys. Rev. D_ 89 (6): 064023.
https://doi.org/10.1103/PhysRevD.89.064023.

Klöckner, Andreas, Nicolas Pinto, Yunsup Lee, B. Catanzaro, Paul Ivanov, and
Ahmed Fasih. 2012. “PyCUDA and PyOpenCL: A Scripting-Based Approach to GPU
Run-Time Code Generation.” _Parallel Computing_ 38 (3): 157–74.
https://doi.org/10.1016/j.parco.2011.09.001.

LIGO Scientific Collaboration. 2018. “LIGO Algorithm Library - LALSuite.” free
software (GPL). https://doi.org/10.7935/GT1W-FZ16.

Prix, Reinhard. 2009. “Gravitational Waves from Spinning Neutron Stars.” In
_Neutron Stars and Pulsars_ , edited by Werner Becker, 357:651–85. Astrophys.
Space Sci. Lib. Berlin Heidelberg: Springer.
https://doi.org/10.1007/978-3-540-76965-1_24.

Riles, Keith. 2017. “Recent searches for continuous gravitational waves.”
_Mod. Phys. Lett. A_ 32 (39): 1730035.
https://doi.org/10.1142/S021773231730035X.

Sieniawska, Magdalena, and Michał Bejger. 2019. “Continuous gravitational
waves from neutron stars: current status and prospects.” _Universe_ 5 (11):
217. https://doi.org/10.3390/universe5110217.

Vousden, W. D., W. M. Farr, and I. Mandel. 2015. “Dynamic temperature
selection for parallel tempering in Markov chain Monte Carlo simulations.”
_Mon. Not. Roy. Astron. Soc._ 455 (2): 1919–37.
https://doi.org/10.1093/mnras/stv2422.

Wette, Karl. 2020. “SWIGLAL: Python and Octave interfaces to the LALSuite
gravitational-wave data analysis libraries.” _SoftwareX_ 12: 100634.
https://doi.org/10.1016/j.softx.2020.100634.