# Metapopulation model of phage therapy of an acute Pseudomonas aeruginosa lung infection

Rogelio A Rodriguez-Gonzalez, Quentin Balacheff, Laurent Debarbieux,

Jacopo Marchi, Joshua S Weitz

► To cite this version:

Rogelio A Rodriguez-Gonzalez, Quentin Balacheff, Laurent Debarbieux, Jacopo Marchi, Joshua S Weitz. Metapopulation model of phage therapy of an acute Pseudomonas aeruginosa lung infection. mSystems, 2024, 9 (10), pp.e0017124. 10.1128/msystems.00171-24. hal-05215189

HAL Id: hal-05215189

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

Submitted on 19 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Metapopulation model of phage therapy of an acute
## *Pseudomonas aeruginosa* lung infection

Rogelio A. Rodriguez-Gonzalez,¹,² Quentin Balacheff,³ Laurent Debarbieux,⁴ Jacopo Marchi,⁵ Joshua S. Weitz⁵,⁶,⁷

**AUTHOR AFFILIATIONS** See affiliation list on p. 18.

**ABSTRACT** Infections caused by multidrug resistant (MDR) pathogenic bacteria are a global health threat. Bacteriophages ("phage") are increasingly used as alternative or last-resort therapeutics to treat patients infected by MDR bacteria. However, the therapeutic outcomes of phage therapy may be limited by the emergence of phage resistance during treatment and/or by physical constraints that impede phage-bacteria interactions in vivo. In this work, we evaluate the role of lung spatial structure on the efficacy of phage therapy for *Pseudomonas aeruginosa* infections. To do so, we developed a spatially structured metapopulation network model based on the geometry of the bronchial tree, including host innate immune responses and the emergence of phage-resistant bacterial mutants. We model the ecological interactions between bacteria, phage, and the host innate immune system at the airway (node) level. The model predicts the synergistic elimination of a *P. aeruginosa* infection due to the combined effects of phage and neutrophils, given the sufficient innate immune activity and efficient phage-induced lysis. The metapopulation model simulations also predict that MDR bacteria are cleared faster at distal nodes of the bronchial tree. Notably, image analysis of lung tissue time series from wild-type and lymphocyte-depleted mice revealed a concordant, statistically significant pattern: infection intensity cleared in the bottom before the top of the lungs. Overall, the combined use of simulations and image analysis of *in vivo* experiments further supports the use of phage therapy for treating acute lung infections caused by *P. aeruginosa*, while highlighting potential limits to therapy in a spatially structured environment given impaired innate immune responses and/or inefficient phage-induced lysis.

**IMPORTANCE** Phage therapy is increasingly employed as a compassionate treatment for severe infections caused by multidrug-resistant (MDR) bacteria. However, the mixed outcomes observed in larger clinical studies highlight a gap in understanding when phage therapy succeeds or fails. Previous research from our team, using *in vivo* experiments and single-compartment mathematical models, demonstrated the synergistic clearance of acute *P. aeruginosa* pneumonia by phage and neutrophils despite the emergence of phage-resistant bacteria. In fact, the lung environment is highly structured, prompting the question of whether immunophage synergy explains the curative treatment of *P. aeruginosa* when incorporating realistic physical connectivity. To address this, we developed a metapopulation network model mimicking the lung branching structure to assess phage therapy efficacy for MDR *P. aeruginosa* pneumonia. The model predicts the synergistic elimination of *P. aeruginosa* by phage and neutrophils but emphasizes potential challenges in spatially structured environments, suggesting that higher innate immune levels may be required for successful bacterial clearance. Model simulations reveal a spatial pattern in pathogen clearance where *P. aeruginosa* are cleared faster at distal nodes of the bronchial tree than in primary nodes. Interestingly, image analysis of infected mice reveals a concordant and statistically significant pattern: infection intensity clears in the bottom before the top of the lungs. The combined use

**Editor** Oleg A. Igoshin, Rice University, Houston, Texas, USA

Address correspondence to Jacopo Marchi, jmarchi@umd.edu, or Joshua S. Weitz, jsweitz@umd.edu.

The authors declare no conflict of interest.

See the funding table on p. 18.

**Received** 6 February 2024

**Accepted** 24 July 2024

**Published** 4 September 2024

Copyright © 2024 Rodriguez-Gonzalez et al. This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license.
of modeling and image analysis supports the application of phage therapy for acute
*P. aeruginosa* pneumonia while emphasizing potential challenges to curative success in
spatially structured *in vivo* environments, including impaired innate immune responses
and reduced phage efficacy.

**KEYWORDS** bacteriophage therapy, microbial ecology, mathematical modeling, infectious disease, lung infection, *Pseudomonas aeruginosa*, virology, antimicrobial agents, antibiotic resistance, innate immunity

The antimicrobial resistance crisis poses a major burden to global health, with ~1.3 million deaths attributed to drug-resistant bacteria in 2019 alone (1). The discovery of new antibiotics has slowed (2), with only two new chemical classes of antibiotics approved since 2017. The antimicrobial resistance crisis has also catalyzed the search for alternative therapies, including phage therapy. Phage therapy conventionally uses bacteriophage (phage) to target and kill specific bacterial pathogens. However, the emergence of phage resistance and/or physical constraints hindering phage propagation can limit phage treatment efficacy.

The ability of phage to kill bacteria is often assessed through *in vitro*, liquid medium assays. In these well-mixed systems, phage can efficiently eliminate both Gram-positive (3) and Gram-negative (3, 4) bacterial pathogens. In contrast, phages are less efficient at clearing bacterial populations with complex, spatial structure (5–8), especially well-established and mature biofilms (5, 7, 8). The spatial heterogeneity of biofilms limits phage propagation and biofilm clearance (9, 10). For instance, slow-growing bacteria living at the center of a biofilm restrict phage propagation and impede biofilm elimination (9). Furthermore, clusters of phage-resistant bacteria can protect susceptible bacteria by blocking phage infection in mixed-strain biofilms (10).

Applications of phage therapy *in vivo* and in the clinic require consideration of the context, including the spatial structure and host immune responses. For example, a three-phage cocktail targeting *Escherichia coli* in the murine gut resulted in the coexistence of phage and bacteria (11). This study also found a significantly higher relative abundance of phage in the luminal part of the ileum than in the mucosal part (11), suggesting that the intestinal mucosa serves as a spatial refuge for bacteria. Further supporting the spatial refuge hypothesis, an independent study found that high T7 phage and *E. coli* titers coexist in the mice colon for up to 3 weeks, with 80% of isolated *E. coli* colonies susceptible to parental phage (12). The host immune responses also play a significant role in shaping phage therapy outcomes *in vivo*. For instance, synergistic interactions between phage and neutrophils enabled clearance of *P. aeruginosa* lung infections in WT mice (13). In contrast, phage therapy failed to eradicate *P. aeruginosa* infection in neutropenic mice (13, 14). Moreover, phage immunogenicity assays indicate that phages do not significantly increase cytokine production and are often well-tolerated by mammalian hosts (13, 15). However, prolonged treatment (~6 months) with intravenously administered phage might induce neutralizing antibody responses against phage (16).

Computational models and theory provide a route to evaluate the impacts of immune
responses on the potential efficacy of phage therapy. For example, influential model-
ing work integrating immune mechanisms into phage therapy (17, 18) assessed the
combined effects of phage life history traits and immune response intensity on bacterial
control. These models assumed that the immune response could grow unconstrained
and that bacteria cannot evade immune responses. As a result, these models predicted
scenarios where the host immune system alone inevitably controls the bacterial infection
(17, 18). Recent work incorporating more nuanced models revealed how phage and a
sufficiently effective innate immune system work synergistically to clear infections (13,
19). Such "immunophage synergy" may also apply to phage-antibiotic combination
therapy (20). It is important to note that these models (13, 17–19) do not explicitly
represent spatial dynamics, a critical component when modeling phage therapy *in vivo*.
Spatially explicit models can be used to assess infection dynamics and therapeutic outcomes in vivo (21). For instance, spatial structure limits phage dispersion and enables spatial refuges for bacteria to survive phage infection, facilitating the coexistence of phage and bacteria (22, 23). In the context of phage-antibiotic combination therapy, an individual-based model (IBM) found that double-resistant mutants emerged when phage dispersion was restricted and antibiotics were heterogeneously distributed. This led to drug-free spatial refuges where bacteria replicate and acquire phage- and antibiotic-resistant mutations (22). Moreover, an IBM of phage treatment of a mixed-strain biofilm identified that fitness costs associated with phage resistance, coupled with spatial protection of susceptible cells vs phage infection by clusters of resistant bacteria, facilitated the long-term coexistence of susceptible and resistant strains (10).

Spatial organization into subpopulations has been proposed as a key factor in fostering the persistence of predator–prey interactions (24). Metapopulation models represent environments as interconnected patches where local subpopulations reside and interact via individuals moving among the patches. A metapopulation model of the host–parasitoid pair, *Callosobruchus chinensis* and *Anisopteromalus clalandrae*, found that the persistence of the host–parasitoid interaction increased with spatial subdivision and with the number of patches (24). In a host plant–pathogen metapopulation model, limited host dispersion led to increased pathogen diversity, while increased host dispersion homogenized the dynamics of local populations and decreased pathogen diversity (25). Additionally, a metapopulation model exploring pathogen circulation in healthcare settings found that pathogen control via host sanitation, rather than environment sanitation, led to more rapid elimination of the circulating pathogen (26).

In this study, we introduce a metapopulation model that integrates lung spatial structure and host immune responses to evaluate their combined impact on phage therapy outcomes for acute pneumonia caused by *P. aeruginosa*. The model simulates ecological interactions between two bacterial strains (phage-susceptible and phage-resistant), phage, and the innate immune response within a network structure resembling the branching pattern of the lungs. Our findings reveal that the clearance of a *P. aeruginosa* infection is contingent upon sufficiently active innate immune states and suitable phage life history traits. Moreover, we note that the spatial model requires higher innate immune response levels to increase the chances of therapeutic success, in contrast to a well-mixed model, which predicts infection resolution with lower innate immune levels, highlighting the spatial structure's role in shaping phage therapy outcomes. Throughout our simulations, we observe a spatial pattern wherein infection clears faster at the network's bottom compared to the top, demonstrating robustness across varied distributions of the bacterial inoculum and phage dose. Lastly, our analysis of in vivo mice lung infection data (13) identifies a concordant and statistically significant pattern: a bottom-to-top spatial transition in pathogen clearance, aligning with the predictions of the metapopulation model.

## Model

The metapopulation network structure is based on the geometry of a symmetrical bronchial tree with a dichotomous branching pattern (Fig. 1a). The nodes of the network represent the airways, and the network links represent the branching points on the bronchial tree. We assume bacteria colonize the airways (nodes) and spread across the network through the network links. We consider the structure of the bronchial tree to be perfectly symmetrical, such that the left and right sides of the tree are identical. Hence, in the absence of stochastic effects, the dynamics at every node of the same tree generation, *g*, will be identical. Therefore, we consider only one node in each generation without losing information on the whole network. The final network topology consists of a chain of connected nodes (Fig. 1a-right), where each node represents one airway in each generation (*g*) of the tree. From the second generation, each level comprises a sister node identical to node *g*, which is topologically connected to *g* and *g* − 1.
<figure><img src="image_2.png" /><figcaption>FIG 1 Schematic of the metapopulation network model of phage therapy of a *P. aeruginosa* infection. The structure of the metapopulation network is based on the geometry of a symmetrical bronchial tree with a dichotomous branching pattern (a). The airways of the bronchial tree (left) are represented by the network nodes (a-middle), while the network links represent the branching points of the bronchial tree. For example, the section that connects the trachea to the left and right main bronchi is considered a branching point. We assume that the bronchial tree is symmetrical such that the left and right parts of the tree are identical and so are their dynamics in a deterministic system. Therefore, we only focus on one side of the tree (the dashed line on the middle network) and reduce the number of nodes in our original network to one node per generation. The final network topology of this reduced model consists of 15 connected nodes that form a chain (right). The degree of the network nodes is 4, except for the first and last nodes, which have a degree of 2. In panel b, we show the ecological interactions between phage-susceptible bacteria (*B*<sub>S</sub>), phage-resistant bacteria (*B*<sub>R</sub>), phage (P), and the host innate immune response (*I*) at the node level (left). Phage infects the susceptible strain, while the resistant strain is targeted only by the innate immune response. The immune response grows in the presence of bacteria and targets both bacterial strains. Neutrophils are recruited to the site of the infection from the pulmonary vasculature, while phage and bacteria transfer between connected nodes to spread across the network (right).</figcaption></figure>

To determine the number of generations of the metapopulation network, we estimate the volume of the network in terms of the ~1 mL total lung capacity (TLC) of mice (27). To do so, we calculate the volume of individual airways of different generations using the anatomical airway information (28). Then, we consider the number of airways per generation (i.e., $2^{g-1}$ airways per generation, g) to estimate the network volume. We estimated that 15 generations yield a network volume of 0.9 mL, which does not exceed the mouse TLC.

Nonlinear interactions between bacteria, phage, and neutrophils govern the population dynamics at the node level, corresponding to a single airway at generation g, Fig. 1b-left. By using a system of nonlinear, ordinary differential equations (equations 1–4), we represent the ecological interactions among bacteria, phage, and host innate immunity occurring at the node level (g). We propose a model with two *P. aeruginosa* strains, one of which is phage-susceptible (*B*<sub>S</sub>) and a second which is phage-resistant
$(B_R)$. The growth of bacteria is density-dependent and is limited by the carrying capacity, $K_C$ can be infected and lysed by phage (P), while $B_R$ resists phage infection. Phage, P, replicates inside bacteria and decay in the environment at a rate $\omega$. Neutrophils target and kill both $B_S$ and $B_R$ bacterial strains. The dynamics are as follows:

$$ \frac{dB_{S,g}}{dt} = r_P B_{S,g} \left(1 - \frac{B_{tot,g}}{K_C}\right) \overbrace{\text{Mutation to } B_R}^{B_S \text{ Growth}} + \mu_2 r_A B_{R,g} \left(1 - \frac{B_{tot,g}}{K_C}\right) - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g-1}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} \left(1 - \frac{V_g - 1}{V_g}\right) \left(1 - \delta_{1,g}\right), \quad (1) $$

$$ \frac{dB_{R,g}}{dt} = r_A B_{R,g} \left(1 - \frac{B_{tot,g}}{K_C}\right) \overbrace{\text{Mutation to } B_S}^{B_R \text{ Growth}} + \mu_1 r_P B_{S,g} \left(1 - \frac{B_{tot,g}}{K_C}\right) + \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} \left(1 - \frac{V_g + 1}{V_g}\right) \left(1 - \delta_{N,g}\right) $$

$$ \frac{dB_{R,g}}{dt} = r_A B_{R,g} \left(1 - \frac{B_{tot,g}}{K_C}\right) \overbrace{\text{Mutation to } B_S}^{B_R \text{ Growth}} + \mu_1 r_P B_{S,g} \left(1 - \frac{B_{tot,g}}{K_C}\right) + \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} \left(1 - \frac{V_g - 1}{V_g}\right) \left(1 - \delta_{1,g}\right), \quad (2) $$

$$ \frac{dB_{R,g}}{dt} = r_A B_{R,g} \left(1 - \frac{B_{tot,g}}{K_C}\right) \overbrace{\text{Mutation to } B_S}^{B_R \text{ Growth}} + \mu_1 r_P B_{S,g} \left(1 - \frac{B_{tot,g}}{K_C}\right) + \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} \left(1 - \frac{V_g - 1}{V_g}\right) \left(1 - \delta_{1,g}\right), \quad (2) $$

$$ \frac{dB_{R,g}}{dt} = r_A B_{R,g} \left(1 - \frac{B_{tot,g}}{K_C}\right) \overbrace{\text{Mutation to } B_S}^{B_R \text{ Growth}} + \mu_1 r_P B_{S,g} \left(1 - \frac{B_{tot,g}}{K_C}\right) + \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} \left(1 - \frac{V_g + 1}{V_g}\right) \left(1 - \delta_{N,g}\right) $$

$$ \frac{dP_g}{dt} = \frac{\text{Viral release}}{\beta B_{S,g} F(P_g)} - \frac{\text{Decay}}{\omega P_g} - \frac{\text{Outflux}}{\frac{8D_P}{l_2 g} P_g} + \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} \left(1 - \frac{V_g - 1}{V_g}\right) \left(1 - \delta_{1,g}\right), \quad (3) $$

$$ \frac{dP_g}{dt} = \frac{\text{Viral release}}{\beta B_{S,g} F(P_g)} - \frac{\text{Decay}}{\omega P_g} - \frac{\text{Outflux}}{\frac{8D_P}{l_2 g} P_g} + \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} - \frac{\overbrace{\text{Immune killing}}^{\text{Lysis}} \overbrace{\text{Outflux}}^{\text{Influx from g-1}} \overbrace{\text{Influx from g}}^{\text{Influx from g}}}{1 + \frac{B_{tot,g}}{K_C}} \left(1 - \frac{V_g + 1}{V_g}\right) \left(1 - \delta_{N,g}\right) $$

Immune stimulation

$$ \frac{dI_g}{dt} = \alpha I_g \left(1 - \frac{I_g}{K_I}\right) \left(\frac{B_{tot,g}}{B_{tot,g} + K_N}\right) \left(1 - \frac{\sum_{g=1}^{15} I_g V_g 2^g - 1}{N_{Lung}}\right), \quad (4) $$

In this model, $B_S$ grows at a maximum rate $r_P$, while $B_R$ grows at a maximum rate $r_A$, and the total bacterial density is represented by $B_{tot} = B_S + B_R$. Phages infect and lyse the phage-susceptible bacteria ($B_S$) at a rate $F(P)$ with a burst size of $\beta$. The host innate immune response ($I$) is activated in the presence of bacteria and grows at a maximum rate $\alpha$, while $K_N$ represents the bacterial density at which the growth of the immune response is half-saturated (29).

We model a saturating immune response by including a local density carrying capacity in the immune growth, $K_I$, as it was shown in previous experimental works that the innate immune response cannot grow unconstrained due to a finite number of immune cells that can be produced and tolerated (30, 31). Moreover, we use the term $N_{Lung}$, representing the observed neutrophil count in the mice lungs (31), to globally constrain neutrophil recruitment across the bronchial network. The innate immune response kills both susceptible and resistant bacteria with a maximum killing rate $\varepsilon$. However, bacteria at high densities can activate mechanisms to evade the immune response, such as biofilm formation (32) and the production of virulence factors (33), resulting in reduced immune killing efficiency. We account for this by introducing the saturation at high bacterial load of the immune killing rate (i.e., we scale the immune killing term by the denominator $1 + B_{tot}/K_D$), in agreement with the experimental measurements of in vivo bacteria burden (34, 35), where $K_D$ is the bacterial density at which the immune killing rate is half its maximum.

These local reaction terms in each node in our metapopulation model build on previous works that established an "immunophage synergy" regime (13, 19) in which both phage and the host innate immune response are required to eliminate the bacterial
infection; i.e., the infection persists when only one or neither is present. The immuno-
phage synergy regime has been previously characterized mathematically (19), where
immune saturation and bacterial immune evasion were identified as key ingredients
leading to such a regime. Immunophage synergy was then observed *in vivo* by infecting
mice with varying immune statuses with pathogenic bacteria and treating them with
phage (13).

Bacteria and phage can hop between connected nodes (Fig. 1b-right). The rate at which they hop is determined by the time, $\tau_g$, required for each species to cross half of the airway length ($l_g$) via diffusion $D$. The hopping rate is calculated as the reciprocal of these times, $\frac{1}{\tau_g}$. In the model, influx and outflux terms govern the transfer of bacteria and phage between neighboring nodes. We assume a homogeneous spreading between airways such that the links between nodes are non-weighted and the flux from node $j$ to node $g$ is inversely proportional to the degree of node $j$ (i.e., $\rho_j = 1/d_j$). Additionally, the number of bacteria and phage transferred from a neighbor node $j$ to a local node $g$ is rescaled to the local density, i.e., $\frac{B_{S,j}V_j}{V_g}$. Note that we multiply by 2 the influx terms from $g+1$ into $g$ since node $g$ is linked to both the main node $g+1$ and its sister node in the reduced network topology (Fig. 1a-right), although we do not need to explicitly model the dynamics in sister nodes by virtue of the branching tree symmetry.

Changes in the concentration of mucin lining the airways impact bacterial motility (36) and phage diffusion (37, 38), thereby affecting their hopping rate throughout the bronchial network. Text S1 B-C explains in detail how to calculate the diffusion constants and the hopping rate of phage and bacteria across different mucin levels. Additionally, the diffusion values of phage and bacteria can be found in Table S1.

For phage infection, we use a model, $F(P)$, that assumes that spatial heterogeneities inside the mouse lungs might limit phage-bacteria encounters. $F(P) = \tilde{\phi}P^{\sigma}$, where $\tilde{\phi}$ is the nonlinear phage adsorption rate and $\sigma < 1$ is the power-law exponent in the phage infection rate. This heterogeneous mixing model (13) has been used previously to recapitulate phage-bacteria dynamics *in vivo*.

The parameter values used in model simulations, largely chosen leveraging the inference work done in (13), are listed in Table 1. Additional details on the model, network structure, calculation of bacteria and phage transfer rates, and statistical analysis are provided in Text S1.

RESULTS

Spatiotemporal dynamics of phage therapy of acute pneumonia caused by P. aeruginosa

We begin by simulating phage therapy treatment of a *P. aeruginosa* infection in an immunocompetent host. The infection starts by inoculating the host with 10⁶ bacterial cells. For phage treatment, we add 10⁷ phage 2 hours after the beginning of the infection. The bacterial inoculum and the phage dose are uniformly distributed among network nodes, ensuring that each node has the same initial bacterial density of 1.11 × 10⁶ CFU/mL and phage density of 1.11 × 10⁷ PFU/mL, considering a network volume of ~0.9 mL. We simulate baseline neutrophil levels by setting an initial immune density of 4.05 × 10⁵ neutrophils/mL in all network nodes (31). Finally, considering that we model phage treatment dynamics in realistic physical compartments with actual volume differences, we show in Fig. 2 the population dynamics (in densities) of bacteria, phage, and the innate immune response at the node level when their numbers exceed 1 CFU (PFU or cell). Densities below this threshold are not shown.

By analyzing the dynamics that emerge from the phage therapy of a *P. aeruginosa* infection (Fig. 2), we observed that the phage-susceptible bacterial population ($B_S$) initially increased and reached its peak density after 13 hours. The phage population ($P$)
<table><caption>TABLE 1 Parameter values of the metapopulation network model. Concentrations are calculated based on the mice lung volume of 0.9 mL</caption><thead><tr><th>Parameters of the metapopulation model</th><th>Value</th><th>Estimated from</th></tr></thead><tbody><tr><td>r<sub>p</sub>, maximum growth rate of phage-susceptible bacteria</td><td>0.75 h<sup>-1</sup></td><td>P. aeruginosa murine pneumonia model (35)</td></tr><tr><td>r<sub>A</sub>, maximum growth rate of phage-resistant bacteria</td><td>0.675 h<sup>-1</sup></td><td>10% tradeoff between resistance against the phage and growth rate (39)</td></tr><tr><td>K<sub>C</sub>, carrying capacity of bacteria</td><td>1.49 &times; 10<sup>9</sup> CFU/mL</td><td>Assuming ~4 times above the typical bacterial density in wild-type mice 24 hours post-infection, estimated from reference 13</td></tr><tr><td>&beta;, burst size of phage</td><td>100</td><td>Estimated from reference 13</td></tr><tr><td>&omega;, decay rate of phage</td><td>0.07 h<sup>-1</sup></td><td>Estimated from reference 13</td></tr><tr><td>&epsilon;, killing rate parameter of the immune response</td><td>5.47 &times; 10<sup>7</sup> mL/(h cell)</td><td>Set such that &epsilon;K<sub>I</sub> gives the maximum granulocyte killing rate (35)</td></tr><tr><td>&alpha;, maximum growth rate of the immune response</td><td>0.07 h<sup>-1</sup></td><td>Fitting of neutrophil recruitment data (31)</td></tr><tr><td>K<sub>I</sub>, maximum capacity of the immune response</td><td>3.59 &times; 10<sup>6</sup> cell/mL</td><td>Fitting of neutrophil recruitment data (31)</td></tr><tr><td>K<sub>D</sub>, bacteria concentration at which the immune response is half as effective</td><td>6.14 &times; 10<sup>6</sup> CFU/mL</td><td>Corresponds to a lethal dose of about 5.5 &times; 10<sup>6</sup> CFU/ lung</td></tr><tr><td>K<sub>N</sub>, bacteria concentration when the immune response growth rate is half its maximum</td><td>10<sup>7</sup> CFU/mL</td><td>In vitro data of TLR5 response to PAK strain (29)</td></tr><tr><td>B<sub>0</sub>, initial bacterial density (when the inoculum is uniformly distributed among network nodes)</td><td>1.1106 &times; 10<sup>6</sup> CFU/mL</td><td>Total inoculum of 10<sup>6</sup> CFU</td></tr><tr><td>P<sub>0</sub>, initial phage dose (when the dose is uniformly distributed among network nodes)</td><td>1.1106 &times; 10<sup>7</sup> PFU/mL</td><td>Total phage dose of 10<sup>7</sup> PFU</td></tr><tr><td>I<sub>0</sub>, initial immune response</td><td>4.048 &times; 10<sup>5</sup> cell/mL</td><td>Fitting of neutrophil recruitment data (31)</td></tr><tr><td>I<sub>0</sub>, initial immune response (neutropenic mice)</td><td>0 cell/mL</td><td>Assuming no primary innate immunity</td></tr><tr><td>N<sub>Lung</sub>, total number of neutrophils in the lungs</td><td>3.24 &times; 10<sup>6</sup> cells</td><td>Fitting of neutrophil recruitment data (31)</td></tr><tr><td>&mu;<sub>1</sub>, probability of emergence of phage-resistant mutants per cellular division</td><td>2.85 &times; 10<sup>-8</sup></td><td>Estimated from experimental measurements (40)</td></tr><tr><td>&mu;<sub>2</sub>, probability of emergence of reversible mutants (i.e., from phage-resistant to phage-susceptible) per cellular division</td><td>2.85 &times; 10<sup>-8</sup></td><td>Approximated to the estimates from reference 40</td></tr><tr><td colspan="3"><strong>Heterogeneous mixing (HM) model parameters</strong></td></tr><tr><td>&phi;&#x0303;, nonlinear phage adsorption rate</td><td>1.686 &times; 10<sup>-7</sup> (mL/PFU)&#x0303;h<sup>-1</sup></td><td>Estimated from reference 13</td></tr><tr><td>&sigma;, power-law exponent in the phage infection rate</td><td>0.6</td><td>Estimated from reference 13</td></tr><tr><td colspan="3"><strong>Anatomical information of airways</strong></td></tr><tr><td>V<sub>g</sub>, airway volume</td><td>4.8 &times; 10<sup>-3</sup> to 2.2 &times; 10<sup>-5</sup> mL</td><td>Calculated as the volume of a cylinder (&pi;r<sup>2</sup>h), using anatomical information of mice airways from reference 28</td></tr><tr><td>l<sub>g</sub>, airway length</td><td>0.76–0.06 cm</td><td>Anatomical information of mice airways from reference 28</td></tr><tr><td>Airway diameter</td><td>0.089–0.02 cm</td><td>Anatomical information of mice airways from reference 28</td></tr></tbody></table>

grew together with the $B_S$ population during the first hours of the simulation. The host innate immune response (I) also increased in the presence of bacteria during the first 10 hours of the simulation. Once the phage population reached high-density levels, i.e., ~10<sup>11</sup> PFU/mL, we observed a reduction in the $B_S$ population across the network. The bacterial elimination rate accelerated as the phage reduced the $B_S$ population density to a level that facilitated the innate immune response to control the infection. Although phage-resistant mutants ($B_R$) dynamics led to transient increase above 10<sup>4</sup>/mL in the top nodes of the network in contrast to bottom nodes (due to threshold effects mentioned
<figure><img src="image_3.png" /><figcaption>**FIG 2** Spatiotemporal dynamics of phage therapy of a *P. aeruginosa* lung infection. Population dynamics at the node level for phage (solid yellow line), phage-susceptible bacteria (solid blue line), phage-resistant bacteria (solid orange line), and the host innate immune response (purple solid line). We set the lower bound of the y-axis for each node panel to the $1/V_i$ density level, such that only densities corresponding to species numbers >1 CFU (PFU or cell) are depicted, where $V_i$ is the volume of node i. Simulation configuration is provided in the main text with parameter values in Table 1 and computational details in Text S1.</figcaption></figure>

in the previous paragraph), the $B_R$ population remained at levels that could be controlled by the host innate immune response. The combined effects of phage and neutrophils led to the clearance of the infection in the network after 32 hours (Fig. 2), similar to observed clearance times of 24–48 hours in phage treatment of *P. aeruginosa* in immunocompetent mice (13). Note that despite similar population dynamics between nodes, the infection clearance did not occur simultaneously in all the nodes. The infection resolved first in the bottom (~28 hours) and later in the top nodes (~33 hours).

## Joint action of phage and host innate immunity in bacterial infections

Next, we study how different treatment scenarios impact the bacterial infection spatiotemporal dynamics across the metapopulation network. We test four treatment scenarios that result from the presence (+) or absence (−) of both phage and innate immunity. As before, we inoculate a total number of $10^6$ bacterial cells that are homogeneously distributed in the network. When phage therapy is applied, we inoculate the host with $10^7$ phage (MOI ≈ 10) 2 hours after the beginning of the infection. If the host is immunocompetent, we set an initial immune density of $4.05 \times 10^5$ neutrophils/mL for all network nodes, totaling $3.64 \times 10^5$ neutrophils in the lungs (31). When the host is immunodeficient, we set the immune density to $I = 0$ neutrophils/mL in the network
nodes. The bacterial inoculum and the phage dose are uniformly distributed in the network such that each node has the same initial bacterial density (1.11 × 10⁶ CFU/mL) and phage density (1.11 × 10⁷ PFU/mL), as in the prior section.

When an immunodeficient host was not treated with phage, bacteria grew unimpe-
ded (Fig. S2a), reaching carrying capacity levels (1.5 × 10⁹ CFU/mL) after 10 hours in all
the nodes (Fig. 3a). On the other hand, an active innate immune response delayed the
initial growth of bacteria due to neutrophils killing the bacteria (Fig. 3b). However, in the
absence of another antimicrobial effect, bacteria continued to grow and overwhelmed
the innate immune response, rendering it insufficient for controlling the infection (Fig.
S2b). The phage treatment of *P. aeruginosa* infection in an immunodeficient host led to
the elimination of the *B*<sub>S</sub> population after 50 hours in most network nodes. However, it
took about 90 hours to eliminate the *B*<sub>S</sub> population in the top nodes (Fig. S2c). During
this time, phage-resistant mutants (*B*<sub>R</sub>) emerged and spread throughout the network,
causing the infection to persist (Fig. 3c). Finally, the model predicted therapeutic success
when host innate immunity complemented phage therapy (Fig. 3d). We found this
outcome to be robust to variations of up to 20% in the critical parameter *K*<sub>D</sub> (Fig. S3a)
and to higher initial bacterial (10<sup>7</sup> CFU) and phage (10<sup>8</sup> PFU) inocula (Fig. S4), which are
conditions consistent with *in vivo* treatment of the *P. aeruginosa* strain PAKlumi with the
phage PAK_P1 (13). Model simulations indicated that bacteria are eliminated from the
network after 32 hours. The elimination of the infection followed a spatial pattern that
started at the bottom and continued to the top of the network (Fig. 3e); this pattern
remained robust to variations in the parameter *K*<sub>D</sub> (Fig. S3b).

The spatial pattern of infection clearance arises in part due to volume differences among network nodes, with the bottom nodes being smaller (and having smaller volumes) than the top nodes. Although our model is based on densities, we define a bacterial extinction threshold based on cell counts in a node. Bacterial elimination at the node level occurs when the number of bacterial cells in a node (airway) falls below the 1 CFU threshold. Therefore, if the bottom nodes have smaller volumes than the top nodes, we expect fewer bacterial cells in the bottom nodes than in the top ones so that the extinction threshold would be reached faster at the bottom of the network. The speed at which the extinction threshold is reached depends on bacteria falling below a critical level ($K_D$), after which the immune system removes bacteria at a fixed rate ($\varepsilon K_I$) until complete elimination within each node. In Text S1 D, we further explain the factors giving rise to the spatial pattern of infection elimination and derive an approximation for the time to infection clearance at the node level.

Next, we sought to characterize how changes in mucin levels—known for shaping
phage and bacteria propagation rates (36–38)—impact the observed clearance pattern.
The spatial clearance pattern was recapitulated for mucin levels varying from 1% to 4%
(Fig. 3f). This suggests infection clearance is robust to changes in conditions that affect
phage and bacterial propagation rates when the host is fully immunocompetent. Overall,
modeling outcomes show that synergistic interactions between phage and neutrophils
may lead to the resolution of the infection in spatially structured environments.

Varying the initial distribution of bacteria and phage in the metapopulation network

Thus far, we have explored the impacts of phage therapy on a *P. aeruginosa* infection model by uniformly distributing the bacterial inoculum and the phage dose in the network. However, the actual distribution of phage and bacteria in the bronchial tree after intranasal inoculation is not necessarily well-constrained. Hence, we try different forms of allocating the bacteria and phage inocula in the metapopulation network. For example, (i) we uniformly distribute the bacterial inoculum such that every node has the same initial density (uniform distribution), (ii) we distribute the inoculum among the first three nodes of the network (top distribution), (iii) we distribute the inoculum among the last twelve nodes of the network (bottom distribution), and (iv) we only
<figure><img src="image_4.png" /><figcaption>FIG 3 Bacterial dynamics under different phage and innate immune treatments. We simulate four treatment scenarios that result from the presence or absence of both phage and the innate immune response. We show the bacterial dynamics across the metapopulation network when the host is immunodeficient untreated (a) or phage-treated (c). Similarly, we show the bacterial dynamics when the host is immunocompetent untreated (b) or phage-treated (d). The heatmaps depict the progression of the bacterial infection across the network; each row represents a network node, g, while the columns indicate the simulation time (hour). The node color represents the bacterial density at a given time. The yellow regions represent high bacterial density, and the white areas represent infection clearance. When the host is immunocompetent and phage-treated, we zoom in and show the infection clearance pattern (e). Effects of varying the mucin level (1%-4%) on the infection clearance time (f). On the heatmaps, row 1 represents Generation 1 and the top of the lungs, while row 15 represents Generation 15 and the bottom of the lungs. The simulation configuration is provided in the main text with parameter values in Table 1 and computational details in Text S1.</figcaption></figure>
inoculate the first node of the network (i.e., intratracheal instillation). We use the same distribution forms as the bacterial inoculum for the phage dose. In total, we test 16 ways of distributing the bacterial inoculum and the phage dose in the network.

When the bacterial inoculum is uniformly distributed, the model predicted bacterial elimination occurred first at the bottom (28 hours) and later at the top nodes (32 to 34 hours) regardless of how the phage dose was allocated in the network (Fig. 4a). When the first three nodes were inoculated with bacteria, it took ~4 hours to colonize the bottom nodes and around 32 hours to eliminate the infection from the network (Fig. S5, 2nd column from the left). In this case, the delay in the bacterial colonization of bottom nodes led to delayed infection dynamics and pathogen elimination in bottom nodes such that clearance times were more similar between bottom (31 hours) and top nodes (32 to 33 hours), regardless of initial phage allocation (Fig. 4b). When the bacterial inoculum was distributed among the last twelve nodes of the network, it took ~37 hours to clear the infection from the network. In this scenario, there was a longer infection clearance time difference between the top (clearance occurred at 37 hours) and the bottom (clearance occurred at 28 hours) of the network, attributed to a delay in bacteria colonizing the top of the network (Fig. 4c). Finally, the inoculation of bacteria in node 1 (trachea) resulted in the colonization of bottom nodes after ~6 hours and elimination of the infection after 36 hours (Fig. S5, last right-column); difference in infection elimination times was observed between the bottom (~34 hours) and the network's top (~37 hours), Fig. 4d. Overall, we observed different bacterial colonization patterns depending on how the bacterial inoculum was initially distributed in the network (Fig. S5). Varying the allocation of the phage dose in the network did not significantly impact bacterial colonization or infection clearance patterns. As in previous simulations (Fig. 2 and 3e), we found a bottom-to-top pattern of infection resolution across different network distributions of phage and bacteria inocula. The pattern is compatible with the bacterial extinction threshold we set in our model based on the number of bacterial cells falling below the 1 CFU threshold and the node volume differences between the bottom and top nodes. For a detailed explanation of the factors contributing to this spatial pattern of infection elimination, refer to Text S1 D.

## Sufficient neutrophil levels and phage adsorption rates are required for effective clearance of *P. aeruginosa* lung infection

Next, we evaluated how intermediate phage efficacy and host innate immunity states impact metapopulation model outcomes, comparing the results with those from the well-mixed case. To model intermediate immune states, we vary the percentage of neutrophil availability in the lungs from 1% to 100%, where 100% availability corresponds to ~3.24 × 10⁶ lung neutrophils in immunocompetent mice (31). To explore intermediate phage efficacy, we vary the phage adsorption rate across a range from 10⁻⁹ to 10⁻⁶ (ml/PFU)⁶h⁻¹. We simulate the phage treatment of a *P. aeruginosa* infection by inoculating a host with 10⁶ bacterial cells and introducing 10⁷ phage 2 hours after the bacterial infection. To assess the robustness of model predictions, we randomize the initial conditions and try 84 different ways of allocating the bacterial inoculum and the phage dose in the network. Then, we calculate the probability of clearing the infection by simulating the different initial conditions, given a specific immune state and phage adsorption rate.

When neutrophil availability is limited (e.g., 10%-35%), the metapopulation model predicted therapeutic success ranging between 25% and 42% for phage adsorption rates ($\tilde{\phi}$) ranging between ~2.5 × 10⁻⁷ and 10⁻⁶ (ml/PFU)⁶h⁻¹ (Fig. 5). In contrast, the well-mixed model, lacking consideration of the spatial structure, predicts infection persistence in the same parameter range. The spatial model predicts therapeutic success ranging between 40% and 100% when $\tilde{\phi} > 10^{-7}(ml/PFU)^6h^{-1}$ and for neutrophil availability >50%. In contrast, infection was always resolved in the well-mixed model within the specified parameter range (Fig. 5, region limited by the solid white line).
<figure><img src="image_5.png" /><figcaption>FIG 4 Infection elimination time given variations in the distribution of the phage dose and bacterial inoculum in the bronchial network. We evaluate different forms of allocating the bacterial inoculum ($B_0$) and the phage dose ($P_0$) among network nodes and calculate the infection clearance time at the node level (g). Distributions include (a) uniform distribution of the bacterial inoculum among all network nodes; (b) distribution of the bacterial inoculum between the first three nodes (c), among the last 12 nodes of the network, or (d) exclusively within the first node of the network. We use the same distribution forms as the bacterial inoculum for the phage dose (different markers inside plots). The simulation configuration is provided in the main text with parameter values in Table 1 and computational details in Text S1.</figcaption></figure>

For example, if we compare the metapopulation model predictions for $\tilde{\phi} = 1.58 \times 10^{-7}$ (ml/PFU)$^\sigma h^{-1}$, neutrophil availability should be >85% to have increased chances of therapeutic success (i.e., >60%), while for the same $\tilde{\phi}$ value, infection always clears in the well-mixed model when neutrophil availability is >55%. These findings highlight the role of spatial structure in shaping phage therapy outcomes and suggest that higher levels of host innate immunity might be required to successfully clear a P. aeruginosa infection when spatial constraints are considered.

Within the metapopulation model outcomes, we noticed a trend where the higher the neutrophil availability, the wider the range of $\tilde{\phi}$ values for which the model predicts increased chances of infection clearance. This suggests that despite spatial structure effects that may affect phage-bacteria encounters and the phage adsorption rate, having a robust immune response facilitates the control of the P. aeruginosa infection. On the other hand, low phage efficacy, e.g., $\tilde{\phi} < 10^{-7}$ (ml/PFU)$^\sigma h^{-1}$, may render phage therapy inefficient even when the host is fully immunocompetent, according to predictions from both the metapopulation and well-mixed models (Fig. 5, black region). We interpret these findings as follows. If phage efficacy is low, phage and bacteria can coexist even in the presence of the immune system. The conditions leading to coexistence depend on the relative balance between phage infection, immune killing, and bacterial growth, as demonstrated in previous theoretical work for a well-mixed model
<figure><img src="image_6.png" /><figcaption>**FIG 5** Probability of therapeutic success given intermediate phage efficacy and host innate immune levels. The probability of therapeutic success is examined given the variation in the percentage of neutrophils available (x-axis, from 1% to 100%) and phage adsorption rate (y-axis, from $10^{-9}$ to $10^{-6}$ (ml/PFU)$^{0}h^{-1}$). Therapeutic success is measured based on 84 different initial distributions of the phage and bacterial inoculum. The colored regions represent a $p > 0$ of clearing the infection, while black regions represent failure to clear the infection, i.e., a $p = 0$ of therapeutic success. The white solid line contours the region of infection clearance predicted by the well-mixed model. The simulation configuration is provided in the main text with parameter values in Table 1 and computational details in Text S1. Here, mucin is set to a 2.5% level.</figcaption></figure>

(19). According to the current metapopulation model, below a phage adsorption rate of $10^{-7}$ (ml/PFU)$^{0}h^{-1}$, the infection cannot be cleared as phage and susceptible bacteria coexist, and eventually, bacteria reach the carrying capacity due to the proliferation of phage-resistant bacteria, even in the presence of a fully active immune response. This coexistence regime is depicted in Fig. S7, where we simulated our metapopulation model for a $\tilde{\phi} = 9 \times 10^{-8}$ (ml/PFU)$^{0}h^{-1}$ just below the infection clearance transition, without allowing for phage resistance. Yet, phage, susceptible bacteria, and the innate immune system enter a coexistence state, and the infection is not cleared. Note that we also explored how variations in mucin and innate immune levels impact phage therapy outcomes. The exploratory analysis confirms the synergistic elimination of the infection by phage and neutrophils given sufficiently active immune states that are weakly impacted by mucin level variation (Fig. S8; Text S1 E).

Analysis of *in vivo* *P. aeruginosa* infection data confirms a bottom-to-top infection clearance pattern

Our finding of a systematic bottom-to-top transition pattern supported by multiple simulations prompted us to examine the spatial patterns of clearance within *in vivo* imaging data from phage treatment of *P. aeruginosa*-infected mice. We analyzed the bioluminescence infection signal from images (Fig. 6) generated in a previous phage therapy study (13). The data set included images of 13 mice, comprising both wild-type (WT) ($N = 4$) and lymphocyte-deficient $Rag2^{-/}-Il2rg^{-/}-$ ($N = 9$) mice groups. Following
<figure><img src="image_7.png" /><figcaption>**FIG 6** *In vivo* *P. aeruginosa* murine pneumonia data. Images depict infection dynamics of *P. aeruginosa* strain PAKlumi *in vivo* during 72 hours. Data are for two mice groups, WT (N = 4; mouse #1 to 4) and *Rag2<sup>-/-</sup> Il2rg<sup>-/-</sup> (N = 9; mouse #5 to 13). The bioluminescence signal represents the intensity of the infection in different mouse regions—a proxy for bacterial densities. A pixel intensity of 1 represents the highest bacterial density, while 0 represents the limit of detection. The white dashed line separates the upper and lower compartments of the mouse respiratory system. The orange and green boxes highlight the approximate time when the total intensity signal drops below the intensity threshold in the lower and upper compartments, respectively. Mice were inoculated with 10<sup>7</sup> *P. aeruginosa* cells, and 2 hours after the bacterial inoculation, mice were treated with phage PAK_P1 (10<sup>8</sup> PFU).</figcaption></figure>

*P. aeruginosa* strain PAKlumi infection (10⁷ CFU), mice were treated with phage PAK_P1 (10⁸ PFU) 2 hours later. We split individual mouse images into two compartments (dashed white line in Fig. 6). The top compartment encompassed the nose and throat,
<figure><img src="image_8.png" /><figcaption>**FIG 7** Time series of the total intensity signal and the infection clearance analysis using *in vivo P. aeruginosa* murine pneumonia data. We show the time series of total intensity signals for the upper and lower compartments of 13 mice (a). The total intensity of one compartment is calculated by adding the pixel intensity values from all pixels making up a compartment. The black dashed line represents the intensity threshold (value of 3) below which the total intensity signal clears. We calculate the time to infection resolution for the upper and lower compartments (b). We use data from 13 mice, including WT (N = 4) and *Rag2<sup>-/-</sup> Il2rg<sup>-/-</sup> (N = 9) mice groups. We used the one-sided Wilcoxon signed rank test to compare the infection clearance time difference between the upper and lower compartments (b).</figcaption></figure>

while the bottom compartment included the lungs. To quantify infection intensity within each compartment, we computed the total intensity signal by summing the pixel intensity values of all pixels within a compartment (Fig. 7a). This total intensity served as a proxy for the bacterial infection status inside the mouse, where higher values indicated higher bacterial density. To determine infection signal clearance in each compartment, we established an intensity threshold of 3, below which we considered the total intensity signal as cleared. In Text S1 F, we include comprehensive details on the image analysis methodology.

Using the time series of imaged mice, we tracked the progression of the infection in the two compartments (Fig. 7a). We observed a more rapid decrease in the intensity signal below the threshold in the bottom compartment (~6–8 hours after the onset of infection) compared to the top compartment (~8–48 hours after the onset of infection). The observation implies a faster clearance of infection in the lower respiratory tract than in the upper respiratory tract of mice. Upon calculating the infection clearance time for both compartments, we found that the intensity signal cleared faster in the bottom compartment than in the top compartment. Specifically, the median time to infection clearance was 7 hours for the bottom compartment and 19 hours for the top compartment, yielding a statistically significant difference with a P-value of 0.0133 (Fig. 7b). This spatial pattern of infection resolution qualitatively aligns with predictions made by the metapopulation model, especially under conditions involving active innate immunity and phage therapy.

DISCUSSION

In this study, we developed a metapopulation model of phage therapy of a *P. aeruginosa* lung infection. We modeled the ecological interactions between phage, two *P. aeruginosa* strains (phage-susceptible and phage-resistant), and the host innate immune response at the airway level. The dynamics arising from the phage therapy of a *P. aeruginosa* infection in an immunocompetent host indicates that phage therapy successfully eliminates bacterial infections when the host innate immune responses complement phage treatment. The therapeutic outcome was robust to heterogeneous
distributions of the phage dose and bacterial inoculum in the network and to variations in mucin levels. However, limited neutrophil availability and a low phage adsorption rate may negatively impact phage treatment efficacy. Moreover, the metapopulation model required higher innate immune response levels to increase the likelihood of therapeutic success, compared to a well-mixed model that achieved infection resolution at lower innate immune levels, highlighting the role of spatial structure in shaping phage therapy outcomes. Across our model simulations, we observed a spatial pattern of infection clearance, with infection resolving faster in the bottom nodes compared to the top nodes. This pattern arises partly due to the dynamics crossing a critical, lower threshold via a volumetric effect influenced by the smaller number of bacterial cells per airway in the bottom vs top nodes. Finally, by analyzing *in vivo* data of mice infected with *P. aeruginosa* and treated with phage, we observed that bacterial infections clear faster in the lower than the upper respiratory tract. This result is qualitatively similar to the bottom-to-top infection elimination pattern predicted by the model. We note that the imaging data for these mice are coarse-grained, preventing us from compartmentalizing the images at an airway-level resolution corresponding to the model scale. Therefore, the mismatch in resolution between the data and the model simulations allows only for a qualitative comparison of the spatiotemporal patterns. Despite this limitation, we found a similar directionality in pathogen elimination in both the data and the model under comparable conditions of immune and phage effects.

The therapeutic outcome recapitulated in our model is consistent with *in vivo* phage therapy results, where phage efficiently controlled respiratory bacterial infections in both humans (41–43) and mice (13, 44–46). Notably, we found that the therapeutic outcome was robust to different distributions of the phage dose across the bronchial network, including heterogeneous and homogeneous dose distributions. This is consistent with outcomes from phage delivery strategies using phage-loaded microparticles to distribute phages throughout the lung, effectively reducing infections caused by *P. aeruginosa* (47) and *S. aureus* (48). Furthermore, Delattre *et al*. developed a computational model that recapitulated phage–bacteria interactions *in vivo* and found that varying the route of administration of the phage dose, i.e., intratracheally or intravenously, did not impact the therapeutic outcomes with both courses producing similar effects (44). While variations in bacterial inoculum allocation impacted colonization patterns within the network, they did not disrupt the therapeutic outcome. These collective findings suggest that therapeutic success might be robust to initial phage and bacterial distributions in the lungs, though quantitative differences can arise with mismatched distributions.

Immunophage synergy does have its limits, whether *in vivo* (13), in mean field models (19), or as shown here, in metapopulation models. For example, a weak innate immune response may negatively impact phage treatment efficacy. Our model predicted a minimum of 45% of neutrophil availability to start clearing bacteria from the lungs, although higher innate immunity levels are needed to increase the likelihood of therapeutic success. This is consistent with *in vivo* experiments showing that phage therapy failed to control a *P. aeruginosa* infection when mice had an impaired innate immune system, e.g., when mice were either MyD88-deficient (13) or neutropenic (13, 14). On the other hand, phage therapy successfully eliminated the infection when mice were immunocompetent (13, 14, 45). Future phage therapy studies and clinical trials should examine the importance of the host immune status in determining treatment efficacy.

The metapopulation network model of phage–bacteria–immune dynamics has caveats. For example, the network structure of the model is based on a symmetrical bronchial tree, while the mouse bronchial tree is asymmetric (49). Hence, future model extensions could consider modeling asymmetric branching structures and their effects on phage and bacteria propagation across the bronchial tree. By doing so, we can evaluate if tree asymmetries facilitate spatial refuges for bacteria to protect themselves from phage attacks. A critical defense mechanism of the respiratory system is mucociliary clearance (MCC) (50), where motile cilia of the airway epithelium transport a mucus
layer out of the lungs, carrying particles trapped in the mucus layer. Currently, our model does not address the (stochastic) mechanism by which bacteria first infect the lungs, the mucosal layer's involvement in this process, or the potential translocation of bacteria between different body sites (51). In the future, we could extend the model to include MCC effects to assess how MCC impacts bacterial colonization, infection clearance patterns, and neutrophil transport within the airways. Our model does not account for the interactions between pathogenic *P. aeruginosa* and the lung microbiome (52). Furthermore, we do not consider the effects of other innate effector cells, such as macrophages, the predominant immune cells in the mice lungs before infection begins (31). The role of alveolar macrophages in the activation of neutrophils, phagocytosis, and the possible interactions these cells have with phage (53, 54) should be considered for future model extensions.

Phage therapy has emerged as an alternative to chemical antibiotics to treat infections by MDR bacteria. However, spatial structure effects of *in vivo* environments where phage and bacteria interact have not been fully addressed. Here, we used a metapopulation network model to provide further support that synergistic interactions between phage and neutrophils can lead to bacterial elimination within the spatially structured lung environment. However, higher innate immune levels could be required to increase the likelihood of therapeutic success in a spatially structured environment compared to a well-mixed system. Our modeling outcomes further demonstrated the importance of the host immune status and phage life history traits, such as phage adsorption rate, in shaping the therapeutic outcomes. In the future, it will be important to obtain high-resolution *in vivo* temporal data and exploit it to achieve fully quantitative model-data comparisons. Such advances will be key to improving the predictive power of *in vivo* phage therapy models at finer scales. The emergence of bottom-to-top spatiotemporal infection clearance patterns also suggests that extending *in vitro* phage therapy models may help guide therapeutic treatment development in other *in vivo* contexts.

# MATERIALS AND METHODS

## Parameter estimation

The parameter values used in the simulations of the metapopulation model are shown in Table 1. Much of the parameter estimation was carried out in previous work (see the "Parameter Estimation" section in reference 13). Airway anatomical information, including airway length and diameter, was obtained from (28). Bacteria and phage diffusion constants can be found in Table S1.

## Model simulation

Model equations (equations 1-4) are numerically integrated using ODE45 in MATLAB R2020b. By doing so, we obtain the population dynamics of phage-susceptible and phage-resistant bacteria, phage, and the innate immune response. We assume bacteria become extinct in a node when the number of bacteria drops below the extinction threshold of 1 CFU. This applies to both $B_S$ and $B_R$. When bacterial counts drop below the extinction threshold in a node, we set the bacterial density to 0 CFU/mL.

# ACKNOWLEDGMENTS

We thank Jeremy Seurat for his insights on the statistical analysis part of the manuscript. We thank Stephen Beckett for his feedback on the metapopulation network structure and the manuscript.

The work was supported by the National Institutes of Health (R01AI146592 to J.S.W., L.D.). J.S.W. was supported, in part, by the Chaires Blaise Pascal program of the Île-de-France region. R.A.R.-G. was supported, in part, by the National Council of Humanities, Science, and Technology of Mexico.
## AUTHOR AFFILIATIONS

¹Interdisciplinary Graduate Program in Quantitative Biosciences, Georgia Institute of Technology, Atlanta, Georgia, USA

²School of Biological Sciences, Georgia Institute of Technology, Atlanta, Georgia, USA

³CHU Félix Guyon, Service des maladies respiratoires, La Réunion, France

⁴Department of Microbiology, Institut Pasteur, Université Paris Cité, CNRS UMR6047, Bacteriophage Bacterium Host, Paris, France

⁵Department of Biology, University of Maryland, College Park, Maryland, USA

⁶Department of Physics, University of Maryland, College Park, Maryland, USA

⁷Institut de Biologie, École Normale Supérieure, Paris, France

## AUTHOR ORCIDs

Rogelio A. Rodriguez-Gonzalez http://orcid.org/0000-0002-5460-4999

Laurent Debarbieux http://orcid.org/0000-0001-6875-5758

Jacopo Marchi http://orcid.org/0000-0003-1192-1203

Joshua S. Weitz http://orcid.org/0000-0002-3433-8312

## FUNDING

<table><thead><tr><td>Funder</td><td>Grant(s)</td><td>Author(s)</td></tr></thead><tbody><tr><td>HHS | National Institutes of Health (NIH)</td><td>R01 AI146592</td><td>Laurent Debarbieux<br>Joshua S. Weitz</td></tr><tr><td>Consejo Nacional de Ciencia y Tecnología (CONACYT)</td><td></td><td>Rogelio A. Rodriguez-Gonzalez</td></tr></tbody></table>

## AUTHOR CONTRIBUTIONS

Rogelio A. Rodriguez-Gonzalez, Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review and editing | Quentin Balacheff, Conceptualization, Investigation | Laurent Debarbieux, Conceptualization, Data curation, Funding acquisition, Methodology, Writing – review and editing | Jacopo Marchi, Conceptualization, Investigation, Methodology, Supervision, Writing – review and editing | Joshua S. Weitz, Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – review and editing

## DATA AVAILABILITY

The code and data used to simulate the metapopulation model, perform the image analysis, and generate the main figures, as well as the supplemental figures, can be found in the GitHub repository at https://github.com/RogerRln/metapop_lung.

## ADDITIONAL FILES

The following material is available online.

## Supplemental Material

**Supplemental material (mSystems00171-24-s0001.pdf).** Supplemental text, table, and figures.

## REFERENCES

1. Murray CJL, Ikuta KS, Sharara F, Swetschinski L, Robles Aguilar G, Gray A, Han C, Bisignano C, Rao P, Wool E, et al. 2022. Global burden of bacterial antimicrobial resistance in 2019: a systematic analysis. The Lancet 399:629–655. https://doi.org/10.1016/S0140-6736(21)02724-0

2. World Health Organization. 2022. 2021 antibacterial agents in clinical and preclinical development: an overview and analysis. World Health Organization.
3. Pallavali RR, Degati VL, Lomada D, Reddy MC, Durbaka VRP. 2017. Isolation and in vitro evaluation of bacteriophages against MDR-bacterial isolates from septic wound infections. PLoS One 12:e0179245. https://doi.org/10.1371/journal.pone.0179245

4. Chibani-Chennoufi S, Sidoti J, Bruttin A, Kutter E, Sarker S, Brüssow H. 2004. *In vitro* and *in vivo* bacteriolytic activities of *Escherichia coli* phages: implications for phage therapy. Antimicrob Agents Chemother 48:2558-2569. https://doi.org/10.1128/AAC.48.7.2558-2569.2004

5. El-Atrees DM, El-Kased RF, Abbas AM, Yassien MA. 2022. Characterization and anti-biofilm activity of bacteriophages against urinary tract enterococcus faecalis isolates. Sci Rep 12:13048. https://doi.org/10.1038/s41598-022-17275-z

6. Rahman M, Kim S, Kim SM, Seol SY, Kim J. 2011. Characterization of induced *Staphylococcus aureus* bacteriophage SAP-26 and its anti-biofilm activity with rifampicin. Biofouling 27:1087-1093. https://doi.org/10.1080/08927014.2011.631169

7. Chaudhry WN, Concepción-Acevedo J, Park T, Andleeb S, Bull JJ, Levin BR. 2017. Synergy and order effects of antibiotics and phages in killing *Pseudomonas aeruginosa* biofilms. PLoS One 12:e0168615. https://doi.org/10.1371/journal.pone.0168615

8. Morris J, Kelly N, Elliott L, Grant A, Wilkinson M, Hazratwala K, McEwen P. 2019. Evaluation of bacteriophage anti-biofilm activity for potential control of orthopedic implant-related infections caused by *Staphylococcus aureus*. Surg Infect (Larchmt) 20:16-24. https://doi.org/10.1089/sur.2018.135

9. Testa S, Berger S, Piccardi P, Oechslin F, Resch G, Mitri S. 2019. Spatial structure affects phage efficacy in infecting dual-strain biofilms of *Pseudomonas aeruginosa*. Commun Biol 2:405. https://doi.org/10.1038/s42003-019-0633-x

10. Simmons EL, Bond MC, Koskella B, Drescher K, Bucci V, Nadell CD. 2020. Biofilm structure promotes coexistence of phage-resistant and phage-susceptible bacteria. mSystems. 5:10-1128. https://doi.org/10.1128/mSystems.00877-19

11. Lourenço M, Chaffringeon L, Lamy-Besnier Q, Pédron T, Campagne P, Eberl C, Bérard M, Stecher B, Debarbieux L, De Sordi L. 2020. The spatial heterogeneity of the gut limits predation and fosters coexistence of bacteria and bacteriophages. Cell Host Microbe 28:390-401. https://doi.org/10.1016/j.chom.2020.06.002

12. Weiss M, Denou E, Bruttin A, Serra-Moreno R, Dillmann M-L, Brüssow H. 2009. *In vivo* replication of T4 and T7 bacteriophages in germ-free mice colonized with *Escherichia coli*. Virology (Auckl) 393:16-23. https://doi.org/10.1016/j.virol.2009.07.020

13. Roach DR, Leung CY, Henry M, Morello E, Singh D, Di Santo JP, Weitz JS, Debarbieux L. 2017. Synergy between the host immune system and bacteriophage is essential for successful phage therapy against an acute respiratory pathogen. Cell Host Microbe 22:38-47. https://doi.org/10.1016/j.chom.2017.06.018

14. Tiwari BR, Kim S, Rahman M, Kim J. 2011. Antibacterial efficacy of lytic *Pseudomonas* bacteriophage in normal and neutropenic mice models. J Microbiol 49:994-999. https://doi.org/10.1007/s12275-011-1512-4

15. Miernikiewicz P, Dąbrowska K, Piotrowicz A, Owczarek B, Wojas-Turek J, Kicielińska J, Rossowska J, Pajtasz-Piasecka E, Hodyra K, Macegoniuk K, Rzewucka K, Kopciuch A, Majka T, Letarov A, Kulikov E, Maciejewski H, Górski A. 2013. T4 phage and its head surface proteins do not stimulate inflammatory mediator production. PLoS One 8:e71036. https://doi.org/10.1371/journal.pone.0071036

16. Dedrick RM, Freeman KG, Nguyen JA, Bahadirli-Talbott A, Smith BE, Wu AE, Ong AS, Lin CT, Ruppel LC, Parrish NM, Hatfull GF, Cohen KA. 2021. Potent antibody-mediated neutralization limits bacteriophage treatment of a pulmonary *Mycobacterium abscessus* infection. Nat Med 27:1357-1361. https://doi.org/10.1038/s41591-021-01403-9

17. Levin BR, Bull JJ. 1996. Phage therapy revisited: the population biology of a bacterial infection and its treatment with bacteriophage and antibiotics. Am Nat 147:881-898. https://doi.org/10.1086/285884

18. Levin BR, Bull JJ. 2004. Population and evolutionary dynamics of phage therapy. Nat Rev Microbiol 2:166-173. https://doi.org/10.1038/nrmicro822

19. Leung CYJ, Weitz JS. 2017. Modeling the synergistic elimination of bacteria by phage and the innate immune system. J Theor Biol 429:241-252. https://doi.org/10.1016/j.jtbi.2017.06.037

20. Rodriguez-Gonzalez RA, Leung CY, Chan BK, Turner PE, Weitz JS. 2020. Quantitative models of phage-antibiotic combination therapy. mSystems 5:e00756-19. https://doi.org/10.1128/mSystems.00756-19

21. Bull JJ, Christensen KA, Scott C, Jack BR, Crandall CJ, Krone SM. 2018. Phage-bacterial dynamics with spatial structure: self organization around phage sinks can promote increased cell densities. Antibiotics (Basel) 7:8. https://doi.org/10.3390/antibiotics7010008

22. Sousa J de, Rocha EPC. 2019. Environmental structure drives resistance to phages and antibiotics during phage therapy and to invading lysogens during colonisation. Sci Rep 9:3149. https://doi.org/10.1038/s41598-019-39773-3

23. Heilmann S, Sneppen K, Krishna S. 2012. Coexistence of phage and bacteria on the boundary of self-organized refuges. Proc Natl Acad Sci USA 109:12828-12833. https://doi.org/10.1073/pnas.1200771109

24. Bonsall MB, French DR, Hassell MP. 2002. Metapopulation structures affect persistence of predator-prey interactions. J Anim Ecol 71:1075-1084. https://doi.org/10.1046/j.1365-2656.2002.00670.x

25. Papaix J, Burdon JJ, Lannou C, Thrall PH. 2014. Evolution of pathogen specialisation in a host metapopulation: joint effects of host and pathogen dispersal. PLoS Comput Biol 10:e1003633. https://doi.org/10.1371/journal.pcbi.1003633

26. Chen S, Lanzas C, Lee C, Zenarosa GL, Arif AA, Dulin M. 2019. Metapopulation model from pathogen's perspective: a versatile framework to quantify pathogen transfer and circulation between environment and hosts. Sci Rep 9:1694.

27. Irvin CG, Bates JHT. 2003. Measuring the lung function in the mouse: the challenge of size. Respir Res 4:1-9. https://doi.org/10.1186/rr199

28. Counter WB, Wang IQ, Farncombe TH, Labiris NR. 2013. Airway and pulmonary vascular measurements using contrast-enhanced micro-CT in rodents. Am J Physiol Lung Cell Mol Physiol 304:L831-43. https://doi.org/10.1152/ajplung.00281.2012

29. Zhang Z, Louboutin J-P, Weiner DJ, Goldberg JB, Wilson JM. 2005. Human airway epithelial cells sense *Pseudomonas aeruginosa* infection via recognition of flagellin by toll-like receptor 5. Infect Immun 73:7151-7160. https://doi.org/10.1128/IAI.73.11.7151-7160.2005

30. Summers C, Rankin SM, Condliffe AM, Singh N, Peters AM, Chilvers ER. 2010. Neutrophil kinetics in health and disease. Trends Immunol 31:318-324. https://doi.org/10.1016/j.it.2010.05.006

31. Reutershan J, Basit A, Galkina EV, Ley K. 2005. Sequential recruitment of neutrophils into lung and bronchoalveolar lavage fluid in LPS-induced acute lung injury. Am J Physiol Lung Cell Mol Physiol 289:L807-15. https://doi.org/10.1152/ajplung.00477.2004

32. Costerton JW, Stewart PS, Greenberg EP. 1999. Bacterial biofilms: a common cause of persistent infections. Science 284:1318-1322. https://doi.org/10.1126/science.284.5418.1318

33. Gellatly SL, Hancock REW. 2013. *Pseudomonas aeruginosa*: new insights into pathogenesis and host defenses. Pathog Dis 67:159-173. https://doi.org/10.1111/2049-632X.12033

34. Drusano GL, Fregeau C, Liu W, Brown DL, Louie A. 2010. Impact of burden on granulocyte clearance of bacteria in a mouse thigh infection model. Antimicrob Agents Chemother 54:4368-4372. https://doi.org/10.1128/AAC.00133-10

35. Drusano GL, VanScoy B, Liu W, Fikes S, Brown D, Louie A. 2011. Saturation of granulocyte kill of *Pseudomonas aeruginosa* in a murine model of pneumonia. Antimicrob Agents Chemother 55:2693-2695. https://doi.org/10.1128/AAC.01687-10

36. Matsui H, Wagner VE, Hill DB, Schwab UE, Rogers TD, Button B, Taylor RM II, Superfine R, Rubinstein M, Iglewski BH, Boucher RC. 2006. A physical linkage between cystic fibrosis airway surface dehydration and *Pseudomonas aeruginosa* biofilms. Proc Natl Acad Sci USA 103:18131-18136. https://doi.org/10.1073/pnas.0606428103

37. Barr JJ, Auro R, Sam-Soon N, Kassegne S, Peters G, Bonilla N, Hatay M, Mourtada S, Bailey B, Youle M, Felts B, Baljon A, Nulton J, Salamon P, Rohwer F. 2015. Subdiffusive motion of bacteriophage in mucosal surfaces increases the frequency of bacterial encounters. Proc Natl Acad Sci USA 112:13675-13680. https://doi.org/10.1073/pnas.1508355112

38. Joiner KL, Baljon A, Barr J, Rohwer F, Luque A. 2019. Impact of bacteria motility in the encounter rates with bacteriophage in mucus. Sci Rep 9:16427. https://doi.org/10.1038/s41598-019-52794-2
39. Bohannan BJM, Kerr B, Jessup CM, Hughes JB, Sandvik G. 2002. Trade-offs and coexistence in microbial microcosms. Antonie Van Leeuwenhoek 81:107–115. https://doi.org/10.1023/a:1020585711378

40. Luria SE, Delbrück M. 1943. Mutations of bacteria from virus sensitivity to virus resistance. Genetics 28:491–511. https://doi.org/10.1093/genetics/28.6.491

41. Maddocks S, Fabijan AP, Ho J, Lin RCY, Ben Zakour NL, Dugan C, Kliman I, Branston S, Morales S, Iredell JR. 2019. Bacteriophage therapy of ventilator-associated pneumonia and empyema caused by *Pseudomonas aeruginosa*. Am J Respir Crit Care Med 200:1179–1181. https://doi.org/10.1164/rccm.201904-0839LE

42. Law N, Logan C, Yung G, Furr C-LL, Lehman SM, Morales S, Rosas F, Gaidamaka A, Bilinsky I, Grint P, Schooley RT, Aslam S. 2019. Successful adjunctive use of bacteriophage therapy for treatment of multidrug-resistant *Pseudomonas aeruginosa* infection in a cystic fibrosis patient. Infection 47:665–668. https://doi.org/10.1007/s15010-019-01319-0

43. Tan X, Chen H, Zhang M, Zhao Y, Jiang Y, Liu X, Huang W, Ma Y. 2021. Clinical experience of personalized phage therapy against carbapenem-resistant *Acinetobacter baumannii* lung infection in a patient with chronic obstructive pulmonary disease. Front Cell Infect Microbiol 11:631585. https://doi.org/10.3389/fcimb.2021.631585

44. Delattre R, Seurat J, Haddad F, Nguyen T-T, Gaborieau B, Kane R, Dufour N, Ricard J-D, Guedj J, Debarbieux L. 2022. Combination of invivo phage therapy data with *in silico* model highlights key parameters for pneumonia treatment efficacy. Cell Rep 39:110825. https://doi.org/10.1016/j.celrep.2022.110825

45. Debarbieux L, Leduc D, Maura D, Morello E, Criscuolo A, Grossi O, Balloy V, Touqui L. 2010. Bacteriophages can treat and prevent Pseudomonas aeruginosa lung infections. J Infect Dis 201:1096–1104. https://doi.org/10.1086/651135

46. Hua Y, Luo T, Yang Y, Dong D, Wang R, Wang Y, Xu M, Guo X, Hu F, He P. 2017. Phage therapy as a promising new treatment for lung infection caused by carbapenem-resistant *Acinetobacter baumannii* in mice. Front Microbiol 8:2659. https://doi.org/10.3389/fmicb.2017.02659

47. Agarwal R, Johnson CT, Imhoff BR, Donlan RM, McCarty NA, García AJ. 2018. Inhaled bacteriophage-loaded polymeric microparticles

ameliorate acute lung infections. Nat Biomed Eng 2:841–849. https://doi.org/10.1038/s41551-018-0263-5

48. Kalelkar PP, Moustafa DA, Riddick M, Goldberg JB, McCarty NA, García AJ. 2022. Bacteriophage-loaded poly(lactic-co-glycolic acid) microparticles mitigate *Staphylococcus aureus* infection and cocultures of *Staphylococcus aureus* and *Pseudomonas aeruginosa*. Adv Healthc Mater 11:e2102539. https://doi.org/10.1002/adhm.202102539

49. Winkler-Heil R, Hofmann W. 2016. Modeling particle deposition in the Balb/c mouse respiratory tract. Inhal Toxicol 28:180–191. https://doi.org/10.3109/08958378.2016.1148801

50. Nawroth JC, van der Does AM, Ryan (Firth) A, Kanso E. 2020. Multiscale mechanics of mucociliary clearance in the lung. Phil Trans R Soc B 375:20190160. https://doi.org/10.1098/rstb.2019.0160

51. Wheatley RM, Caballero JD, van der Schalk TE, De Winter FHR, Shaw LP, Kapel N, Recanatini C, Timbermont L, Kluytmans J, Esser M, Lacoma A, Prat-Aymerich C, Oliver A, Kumar-Singh S, Malhotra-Kumar S, Craig MacLean R. 2022. Gut to lung translocation and antibiotic mediated selection shape the dynamics of *Pseudomonas aeruginosa* in an ICU patient. Nat Commun 13:6523. https://doi.org/10.1038/s41467-022-34101-2

52. Dickson RP, Erb-Downward JR, Falkowski NR, Hunter EM, Ashley SL, Huffnagle GB. 2018. The lung microbiota of healthy mice are highly variable, cluster by environment, and reflect variation in baseline lung innate immunity. Am J Respir Crit Care Med 198:497–508. https://doi.org/10.1164/rccm.201711-2180OC

53. Jończyk-Matysiak E, Weber-Dąbrowska B, Owczarek B, Międzybrodzki R, Łusiak-Szelachowska M, Łodej N, Górski A. 2017. Phage-phagocyte interactions and their implications for phage application as therapeutics. Viruses 9:150. https://doi.org/10.3390/v9060150

54. Zborowsky S, Seurat J, Balacheff Q, Ecomard S, Nguyen Ngoc Minh C, Titécat M, Evrard E, Rodriguez-Gonzalez RA, Marchi J, Weitz JS, Debarbieux L. 2024. Macrophage-induced reduction of bacteriophage density limits the efficacy of *in vivo* pulmonary phage therapy. bioRxiv:2024.01.16.575879. https://doi.org/10.1101/2024.01.16.575879