Transl Clin Pharmacol. 2023 Dec;31(4):167-190. English.
Published online Nov 22, 2023.
Copyright © 2023 Translational and Clinical Pharmacology
Review

Phage-host-immune system dynamics in bacteriophage therapy: basic principles and mathematical models

Dongwoo Chae
    • Department of Pharmacology, Yonsei University College of Medicine, Seoul 03722, Korea.
Received September 08, 2023; Revised October 17, 2023; Accepted October 19, 2023.

It is identical to the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0/).

Abstract

Phage therapy is progressively being recognized as a viable alternative to conventional antibiotic treatments, particularly in the context of multi-drug resistant bacterial challenges. However, the intricacies of the pharmacokinetics and pharmacodynamics (PKPD) pertaining to phages remain inadequately elucidated. A salient characteristic of phage PKPD is the inherent ability of phages to undergo replication. In this review, I proffer mathematical models that delineate the intricate dynamics encompassing the phage, the host organism, and the immune system. Fundamental tenets associated with proliferative and inundation thresholds are explored, and distinctions between active and passive therapies are accentuated. Furthermore, I present models that aim to illuminate the multifaceted interactions amongst diverse phage strains and bacterial subpopulations, each possessing distinct sensitivities to phages. The synergistic relationship between phages and the immune system is critically examined, demonstrating how the host’s immunological function can influence the requisite phage dose for an optimal therapeutic outcome. A profound understanding of the presented modeling methodologies is paramount for researchers in the realms of clinical pharmacology and PKPD modeling interested in phage therapy. Such insights facilitate a more nuanced interpretation of dose-response relationships, enable the selection of potent phages, and aid in the optimization of phage cocktails.

Keywords
Phage Therapy; Mathematical Model; Population Dynamics

INTRODUCTION

Phage therapy - using bacteriophages (viruses that target bacteria) to counter bacterial infections - has seen a revitalized interest lately [1]. This uptick in attention arises chiefly from the escalating antibiotic resistance dilemma, emphasizing the need for innovative therapeutic approaches [2]. Bacteriophages, often abbreviated to phages, are highly specific, self-amplifying, and adaptable, showing promise in fighting bacterial infections.

Found ubiquitously in nature, there is an estimated 1031 phages on Earth, outnumbering bacteria by tenfold [2]. Their remarkable specificity often limits their target to a particular bacterial species or strain. Such specificity has a double-edged effect: it ensures phage therapy scarcely impacts the host's microbiota, preserving beneficial bacteria. However, it also necessitates exact bacterial pathogen identification and the creation of tailored phage cocktails for effective therapy [3].

When it comes to drug development, evaluating the pharmacokinetic-pharmacodynamic (PKPD) properties is paramount. But for phages, their inherent ability to reproduce makes their PKPD assessment distinct from antibiotics. Mathematical models have emerged as invaluable in studying phage-bacteria dynamics, paving the way for predicting phage therapy outcomes. Leveraging equations that delineate the involved biological processes, these models delve deep into the intricate interactions amongst phages, bacteria, and the host immune defense.

Previous works have been rooted in in vitro time-kill studies [4], animal in vivo PKPD investigations [5], or a combination thereof [6]. Both the synergistic interactions of phage with the immune system [5, 7] and with antibiotics [8] have been rigorously scrutinized through quantitative mathematical frameworks. Furthermore, both web-based [9] and standalone software applications [10] have been devised to predict the intricate dynamics of phage infections. Theoretical explorations have been meticulously conducted to address the co-evolutionary dynamics intertwining bacteria and bacteriophage [11], strategies for optimizing phage cocktails [12], and the emergence of phenotypic resistance [13], among other pertinent subjects. Current endeavors in mathematical modeling persistently illuminate myriad dimensions inherent to phage therapy. This realm undoubtedly represents an avant-garde and burgeoning field of pharmacometrics and systems pharmacology.

Despite ongoing research, no human PKPD studies have been published to date. Currently, 11 trials have been completed, with some findings already published [14-16], and there are six active clinical trials investigating various facets of phage therapy, as recorded on ClinicalTrials.gov. Although there are no FDA-approved phage therapy products, the leading product is in Phase 2 of clinical development. Nevertheless, there are numerous clinical case reports on phage therapy outcomes under emergency or treatment IND conditions [17-20]. The outcomes thus far have been varied, and the safety and efficacy of phage therapy in humans have not yet been conclusively established. Consequently, significant efforts are required to transition phage therapy into routine clinical practice.

In this review, I will introduce the basic principles of phage-bacteria interactions, spotlight the important concepts associated with proliferative and inundative thresholds, distinguish between active and passive phage therapies, and present mathematical models that elucidate phage-host-immune system dynamics. I will then demonstrate how these models can be harnessed to elucidate dose-response relationships and fine-tune phage dosages.

BASIC PRINCIPLES OF PHAGE-BACTERIA INTERACTIONS

For bacteriophages (phages) to kill bacterial cells, they must first adhere to the host bacteria, injecting their genetic content into the cell. From a pharmacokinetic (PK) standpoint, phages are expended after this adsorption. Consequently, the phage loss rate correlates with the bacterial kill rate. This parallels the concept of target-mediated drug disposition (TMDD), where target engagement precedes receptor-mediated endocytosis. Nevertheless, phages have an added ability: replication. There are two core phage replication mechanisms: the lytic and lysogenic cycles [21].

During the lytic cycle, a phage commandeers the host cell's functions to duplicate its DNA and create new virus parts. Once the cell is saturated with phage offspring, it ruptures, dispersing the newborn phage particles to infect other bacteria. This lytic action, which leads to the bacterial host's destruction, is pivotal for therapeutic phages.

Conversely, the lysogenic cycle involves the phage DNA integrating into the bacterial DNA, generating a prophage. This prophage remains latent and duplicates with the host DNA. Under specific triggers, like DNA damage, the prophage may revert to the lytic cycle. While lysogenic phages occasionally have therapeutic applications, they are generally avoided due to risks, such as transferring antibiotic resistance genes to other bacteria. For this discourse, we will primarily concentrate on lytic phages.

MODELING PHAGE INFECTION DYNAMICS

The infection commences with adsorption, where phages attach to receptors on the bacterial cell. The adsorption rate hinges on phage and bacteria concentrations and the phage’s receptor affinity. The constant representing this rate, denoted as k in this article, is vital in mathematical phage-bacteria interaction models.

Following successful adsorption and genome insertion, phage replication ensues, culminating in bacterial cell rupture and progeny release. The term “burst size,” henceforth denoted as b, signifies the number of phages released post-infection, while “latent period” indicates the duration from infection to cell rupture.

Several methods exist to model phage-bacteria dynamics. The predominant technique employs ordinary differential equations (ODEs). However, ODE-based models carry assumptions that can affect their validity:

  • 1. Homogeneous Distribution: Such models often propose a uniform spread of phages and bacteria, a scenario not always reflective of intricate biological systems.

  • 2. Static Parameters: Most models lean on constant values like adsorption rates (k) or burst sizes (b). However, these might fluctuate due to bacterial mutations or shifting environmental conditions.

  • 3. Continuum Assumption: Similar to numerous pharmacokinetic-pharmacodynamic models, these models presuppose that phage and bacteria concentrations are continuous. But, phages and bacteria are distinct units, making these models potentially flawed when phage and bacteria counts are minimal.

Presuming these assumptions are accurate, a rudimentary model delineating phage-host dynamics can be outlined as follows:

(Model I)

dBdt=r·B·1-BBmax-k·P·B(1) 
dPdt=b-1·k·P·B-w·P(2) 

where B represents the bacterial density, P the phage density, r the bacterial growth rate, k the adsorption rate constant, b the burst size, Bmax the carrying capacity, and w the phage decay rate.

To account for the latent period between phage infection and host cell lysis, Model I can be modified into a system of delay differential equations (DDEs).

(Model II)

dSdt=r·S·1-BBmax-k·P·S(3) 
dIdt=k·P·S-k·PL·SL(4) 
dPdt=b·k·PL·SL-k·P·B-w·P(5) 
B=S+I(6) 

where S and I represent sensitive (i.e., uninfected) and infected bacteria. PL and SL denote P and S that existed L time units in the past, with L being the length of the latent period.

Alternatively, transit compartments can be used to approximate the time delay.

(Model III)

dI1dt=k·P·S-ktr·I1(7) 
...
dIndt=ktr·In-1-ktr·In(8) 
dPdt=b·ktr·In-k·P·B-w·P(9) 
B=S+I1+…+In(10) 

For practical purposes, Model I is often sufficient since latent periods of phages used in phage therapy seldom exceeds 1 hour.

UNDERSTANDING PROLIFERATIVE AND INUNDATION THRESHOLDS

Grasping phage PKPD hinges on the notions of proliferative and inundation thresholds [22]. The proliferative threshold is the critical bacterial concentration necessary for a net phage density increase. Conversely, the inundation threshold is the essential phage concentration for a net bacterial density decrease.

To delve into these concepts quantitatively, we refer to a streamlined version of Model I, which overlooks the logistic growth constraint. This model is applicable primarily when B << Bmax.

dBdt=r·B-k·P·B(1’) 

For a net decrease in bacterial concentration, dBdt < 0 must be true. Assuming B > 0, this equates to r > k · P, or P < rk. Therefore, rk is the inundation threshold.

Likewise, the condition for phage replication is gleaned by ensuring dPdt > 0 in Eq. (2). Given P > 0, rearranging Eq. (2) leads to (b − 1) · k · B > w. This makes wb-1·k the proliferation threshold essential for phage propagation.

Fig. 1 illustrates the typical temporal trajectories of lytic phages and their host bacteria. Initially, the bacterial concentration rises as phage density lingers below the inundation threshold. When bacterial density hits the proliferative threshold, there is a sharp uptick in phage density. Over time, as phage density eclipses the inundation threshold, bacterial numbers dwindle. However, the depletion of host bacteria triggers a dip in phage concentration beneath the inundation threshold, reinitiating the cycle.

Figure 1
A schematic diagram illustrating the infection cycle of a typical lytic phage. The bottom graph depicts the expected temporal trajectories of phages and host bacteria, generated using Model I. Parameters used for the simulation are: B(0) = 104 CFU mL−1,P(0) = 104 PFU mL−1,r = 1 h−1,k = 10−7 mL CFU−1 h−1,b = 100, w =1 h−1, and Bmax = 1010 CFU mL−1.

The inundation threshold in phage therapy bears conceptual resemblance to the minimum inhibitory concentration (MIC) observed with antibiotics. However, in phage therapy, even if the initial phage dose does not surpass the inundation threshold, it can still be effective provided there is an adequate bacterial density to foster phage proliferation. This proliferation can eventually elevate the phage density beyond the inundation threshold.

This dynamic sets phage therapy apart from traditional antibiotics: the initial dose is not the sole determinant of the steady-state phage concentration. Instead, phage therapy emphasizes the replicative potential of phages in conjunction with host bacterial density, both of which play pivotal roles in dictating the final phage concentration.

EXPLORATORY SIMULATIONS

Dose-response relationship

Using Model I as a foundation, we explore the dose-response relationship in phage therapy. Fig. 2A and B illustrate that a higher phage dose corresponds to an earlier inhibition of bacterial growth – a predictable outcome. However, what is surprising is that for doses less than or equal to 107 PFU/mL, a higher dose results in a diminished peak phage concentration (Fig. 2A) and an elevated minimum bacterial density (Fig. 2B).

Figure 2
Simulations of dose-response based on Model I. Here, P and B represent phage and bacteria, respectively. Parameters used for the simulations were: B(0) = 104 CFU mL−1,r = 1 h−1,k = 10−7 mL CFU−1 h−1,b = 100, w =1 h−1, and Bmax = 1010 CFU mL−1.

For doses surpassing 107 PFU mL−1, the relationship shifts. Higher doses align with a superior peak phage concentration and a subsequent reduction in bacterial density by t = 6 (Fig. 2B). The U-shaped dose response is a result of rapid bacterial growth inhibition combined with reduced phage multiplication at elevated doses. An essential threshold emerges when the phage dose exceeds a particular point. As depicted in Figure 2B, starting with a phage concentration over 107 PFU mL−1 initiates a dose-dependent bacterial reduction. Therefore, with adequately high doses, the outcome is less reliant on phage multiplication, making the initial dose the primary determinant of final bacterial count.

This pivotal finding suggests 2 distinct therapeutic approaches with phages: active and passive therapy [4]. Active therapy capitalizes on phage multiplication to surpass the inundation threshold, while passive therapy provides a phage dose that already exceeds the threshold. The chosen approach significantly influences the dosing strategy, and the optimal dosage must strike a balance between immediate efficacy and long-term outcomes.

The role of host density

An intriguing facet of phage therapy is its efficacy’s dependence on host density. Elevated host density corresponds with a more substantial peak phage concentration, facilitating more efficient bacterial elimination (Fig. 3A). Theoretically, this means phage therapy is more potent during severe infections, and more effective at infection sites with denser bacterial populations.

Figure 3
The effects of (A) initial host density B(0) and (B) phage elimination rate w on the dynamics of P and B based on Model I. Default parameters used for the sensitivity analysis were: B(0) = 106 CFU mL−1,P(0) = 106 PFU mL−1,r = 1 h−1,k = 10−7 mL CFU−1 h−1,b = 100, w =1 h−1, and Bmax = 1010 CFU mL−1.

Influence of phage pharmacokinetics

A shorter phage half-life results in an elevated minimum bacterial density and a reduced span of the predator-prey cycle (Fig. 3B). Given this, longer-lived phages are more desirable. It is essential to note, however, that a phage’s half-life in its target area might differ from its plasma half-life. Therefore, during preclinical development, it is crucial to measure phage concentrations from the targeted site to acquire accurate data on the phage’s half-life and its precise location of action.

Influence of phage characteristics

Two primary phage parameters – the adsorption rate constant (k) and burst size (b) – play pivotal roles in determining treatment efficacy. Enhancing either k or b leads to a decrease in peak bacterial density. Yet, their implications vary. For simplification, when we exclude phage elimination, simulations indicate a higher k, owing to more rapid bacterial elimination, results in a reduced equilibrium phage titer (Fig. 4A). In contrast, an elevated b produces a heightened equilibrium phage titer (Fig. 4B). Consequently, the rate of bacterial reduction, influenced by the phage titer, diminishes with an increase in k and amplifies with a boost in b.

Figure 4
The effects of (A) adsorption rate constant (k) and (B) burst size (b) on the dynamics of P and B based on Model I. Default parameters used for the sensitivity analysis were: B(0) = 106 CFU mL−1,P(0) = 106 PFU mL−1,r = 1 h−1,k = 10−7 mL CFU−1 h−1,b = 100, w =0 h−1, and Bmax = 1010 CFU mL−1.

This data suggests that when sourcing phages, priority might be given to those with a more considerable burst size (b). Nonetheless, an expanded burst size invariably implies a prolonged latent period [23].

MODELS FOR MULTI-BACTERIA INTERACTIONS

Our assumptions thus far have centered around a uniform bacterial population susceptible to phage. In real-world scenarios, there always exists a subset of bacteria resistant to the phage. To reflect this complexity, bacteria B can be categorized into sensitive (S) and resistant (R) strains:

(Model IV)

dSdt=rS·S·1-cSS·S-cRS·R-k·P·S(11) 
dRdt=rR·R·1-cSR·S-cRR·R(12) 
B=S+R(13) 

where rS and rR depict the growth rate constants of S and R, respectively, and cxy represents the competition coefficient that delineates x’s inhibitory effect on y.

A captivating dynamic emerges, especially when cSR is significantly high (Fig. 5A). An elevated phage dose predominantly eliminates S, prompting S’s reduced inhibition on R. This shift means that higher phage dosages can accelerate the rise of resistance. Such potential outcomes mandate careful consideration during phage dose optimization.

Figure 5
(A) Simulation of Model IV assuming cSS = cRS = cRR = 10−8 CFU−1 mL and cSR = 5 × 10−6 CFU−1 mL. The initial ratio of S:R was set to 99:1 and initial bacterial density to 104 CFU mL−1. Other parameters were set to rS = 1, rR = 2, k = 10−7 mL CFU−1 h−1,b = 100, and w = 0 h−1. (B) Simulation of Model V assuming cnm = 10−8 CFU−1 mL (for n, m = 1, 2, 3), initial B1:B2:B3 = 3,000:3:1, k1 = 10−7 mL CFU−1 h−1,k1 = 10−9 mL CFU−1 h−1,k3 = 10−15 mL CFU−1 h−1,b = 100, and w = 0 h−1.

Expanding on Model IV, an extension can accommodate n subpopulations, each displaying varied phage sensitivities:

(Model V)

dB1dt=r1·B1·1-cn1·Bn-k1·P·B1(14) 
dBndt=rn·Bn·1-cnn·Bn-kn·P·Bn(15) 
B=Bn(16) 

Fig. 5B illustrates the simulation results when n equals 3 subpopulations. Intriguingly, the number of peak bacterial concentrations matches the number of subpopulations incorporated into the model. This characteristic provides information regarding the number of different subpopulations when constructing phage therapy PKPD models.

MULTI-PHAGE INTERACTION MODELS

In multi-phage interaction models, two or more phage strains target the same bacterial strain. These models can help researchers understand how phage strains interact and compete for limited bacterial resources, as well as how phage diversity can impact therapy outcomes. The general form of the multi-phage interaction model is:

(Model VI)

dBdt=r·B·1-BBmax-ki·Pi·B(17) 
dPidt=bi-1·ki·Pi·B-wi·Pi(18) 

where i represents each phage strain, Pi its density, ki the adsorption rate constant, and wi the decay rate of phage strain i.

MULTI-PHAGE, MULTI-BACTERIA INTERACTION MODELS

Model VI can be combined with Model V to yield n subpopulations that differ in sensitivity to the different phages. Such models could be used when modeling the PKPD of phage cocktails.

(Model VII)

dBkdt=rk·Bk·1-BkBmax-kik·Pi·Bk(19) 
dPidt=kbki-1·kik·Pi·Bk-wi·Pi(20) 

where k represents each bacterial subpopulation, Bk its concentration, kik the adsorption rate constant, bki the burst size of phage strain i on bacterial subpopulation k.

The afore-mentioned models have numerous applications in phage therapy research, including:

  • 1) Designing phage cocktails: By simulating the dynamics of multiple phage strains targeting different bacterial strains, these models can help identify the most effective combinations of phages for treating complex bacterial infections.

  • 2) Predicting phage therapy outcomes: Understanding how multiple phages and bacteria interact can improve predictions of therapy outcomes and inform treatment strategies.

  • 3) Investigating coevolutionary dynamics: These models can be used to study the coevolution of phages and bacteria, shedding light on the emergence of bacterial resistance and the potential for phage adaptation in response to resistant bacterial strains.

  • 4) Evaluating the impact of phage therapy on microbiota: By simulating the interactions between multiple phages and bacteria, researchers can assess the potential impact of phage therapy on the overall microbial community and develop strategies to minimize unintended consequences on the host microbiota.

THE HOST IMMUNE SYSTEM AND PHAGE THERAPY

The host immune system is a complex network of cells, tissues, and signaling molecules that work together to protect the body against invading pathogens, including bacteria and viruses. In the context of phage therapy, the immune system can influence the outcome of treatment in several ways:

  • 1) Direct interactions with phages: The immune system can recognize and neutralize phage particles, potentially reducing their efficacy in clearing bacterial infections.

  • 2) Modulation of bacterial clearance: The immune system contributes to bacterial clearance by phagocytosis, the production of antimicrobial peptides, and the activation of other immune cells. The effectiveness of phage therapy can be enhanced or diminished by the host immune response against bacteria.

  • 3) Inflammation and tissue damage: The immune response to bacterial infections can lead to inflammation and tissue damage, which may impact phage therapy outcomes.

To incorporate the host immune system into phage therapy models, researchers can add additional equations and parameters that represent various aspects of the immune response, such as:

  • 1) Phage neutralization: An additional term can be added to the phage dynamics equation to account for the immune system’s capacity to neutralize phage particles. This term may be dependent on the density of phages, as well as the immune cell population.

  • 2) Immune cell dynamics: Equations describing the dynamics of immune cell populations, such as macrophages or neutrophils, can be included in the model. These equations may involve factors such as immune cell recruitment, activation, and proliferation.

  • 3) Interactions between immune cells, phages, and bacteria: The model can be extended to include interactions between immune cells, phages, and bacteria, such as phagocytosis of bacteria by immune cells or the impact of immune cell activation on phage efficacy.

Several mathematical models have been developed that incorporate the host immune system in phage therapy simulations, including:

  • 1) Models investigating the impact of phage neutralization on therapy outcomes: These models can help researchers understand how the host immune system’s ability to recognize and neutralize phage particles can influence the success of phage therapy.

  • 2) Models exploring the synergistic effects of phages and immune cells in bacterial clearance: By simulating the combined action of phages and immune cells against bacteria, these models can provide insights into the optimal conditions for successful phage therapy.

  • 3) Models studying the role of inflammation and tissue damage in phage therapy: These models can help identify strategies to minimize the potential negative effects of immune responses on phage therapy outcomes.

The basic model with phage neutralization

In this model, we explore the neutralization of phage particles by the host’s immune system. We introduce a new variable, I, symbolizing the density of immune cells (like antibodies or phagocytes) that can neutralize phages. The term α · I · P is added to denote the rate of phage neutralization:

(Model VIII)

dBdt=r·B·1-BBmax-k·P·B(21) 
dPdt=b-1·k·P·B-w·P-α·I·P(22) 
dIdt=fB,P,I-δ·I(23) 

Here, α is the neutralization rate constant, while f (B, P, I) describes how immune cell production shifts in response to bacteria (B) and phages (P). δ signifies the decay rate of immune cells.

In the following model, we address the behavior of immune cells, such as neutrophils and macrophages, when faced with bacterial infections. An extra term, g(B, I), embodies the bacterial clearance rate due to immune cells:

(Model IX)

dBdt=r·B·1-BBmax-k·P·B-gB,I(24) 
dPdt=b-1·k·P·B-w·P(25) 
dIdt=hB,P,I-δ·I(26) 

Here, g(B, I) denotes the interaction between bacteria and immune cells, such as phagocytosis. Meanwhile, h(B, P, I) encapsulates the recruitment, activation, and multiplication of immune cells reacting to bacteria (B) and phages (P). δ remains as the decay rate of immune cells.

Diverse dynamic patterns can emerge based on the functional forms assumed for g(B, I) and h(B, P, I). Leung and Weitz [7], under the premise of saturable immune cell killing and immune saturation and assuming δ to be negligible, set the functions as:

gB,I=ε·I·B/(1+BKD)(27) 
hB,P,I=α·I·(1-IKI)·BB+KN(28) 

In the above, ε reflects the killing rate of innate immune response, KD the bacterial concentration at which innate immune response is half saturated, α the maximum growth rate of innate immune response, KI the maximum capacity of innate immune response, and KN the bacterial concentration when innate immune response growth rate is half its maximum.

With this setup, three theoretical outcomes were identified: 1) Joint action of phage and the immune system eradicates bacteria, 2) Phages and bacteria achieve a stable coexistence, 3) Phages get eliminated, leading to the maximal growth of bacteria.

Here, we deal with the simplified version of the model by assuming a steady state level I for the immune system. Saturable bacterial elimination by the immune system is also assumed, in line with Leung and Weitz:

(Model X)

dBdt=r·B·1-BBmax-k·P·B-ε·I·B/(1+BKD)(29) 
dPdt=b-1·k·P·B-w·P(30) 

Without phages, the B equation simplifies to an equation with a constant ε′ (defined as ε · I).

dBdt=r·B·1-BBmax-ε'·B/(1+BKD)29’ 

By setting the change in B over time to zero and comparing the two growth and death functions for B,fB=r·B·1-BBmax and gB=ε'·B/1+BKD, respectively, several insights emerge about bacterial population dynamics and the roles of the immune system and phage.

The function f(B) is an inverted parabola intersecting the abscissa at B = 0 and B = Bmax. The function g(B), on the other hand, is a Michaelian curve approaching the asymptotic maximum of ε · KD as B → ∞. At small B, f(B) ≍ r · B and g(B) ≍ ε′ · B. Hence, when r > ε′, the two functions always intersect at 2 points with B invariably converging to the higher, nonzero solution (Fig. 6A). On the other hand, with ε′r, the 2 functions either intersect at a single point (Fig. 6B) or 3 points (Fig. 6C).

Figure 6
Blue and red curves correspond to f(B) to g(B), respectively. Green dotted curves are series of g(B) under parameter perturbation. (A) The effect of decreasing ε′ from 1 to 0. (B) The effect of increasing ε′ from 1 to 2. (C) The effect of decreasing KD from 1 to 0.1 given ε′ = 2. (D) The bifurcation diagram showing the fixed point(s) of B as KD changes within the range of 0 to 1. Solid and dottetd curves represent stable and unstable fixed points, respectively.

Biological interpretation suggests that under insufficient immune cell kill characterized by ε′ < r, bacteria grow until stabilizing at a certain steady state level BSS with BSS < Bmax. On the other hand, with ε′r, death rate always dominates growth rate at small B. In this case, the magnitude of KD determines whether f(B) and g(B) intersect at a single point or three points.

In essence, the dynamics can exhibit “bi-stability,” meaning the steady state level depends on the system’s initial state. The system can also undergo a “saddle-node bifurcation” at a critical value of KD (Fig. 6D). When KD is above the threshold, zero is the only fixed point of B and is globally stable. Otherwise, 2 additional fixed points emerge, with a lower unstable fixed point (dotted line in Fig. 6D) and a higher stable fixed point (solid line in Fig. 6D). In this case, the initial bacterial density determines the steady state of the system. If the initial bacterial density is lower than the unstable fixed point, B converges to the lower fixed point of zero. However, if the initial bacterial density is higher that the unstable fixed point, bacteria grow and converges to the higher fixed point.

The findings suggest that minor infections with low bacterial density are typically managed by the immune system, while severe infections are not. This matches clinical observations that mild infections often resolve on their own.

Reintroducing the phage, it is assumed that phage is administered at the beginning (t = 0), and the immune system activation has not started. Phages reduce bacterial density. At a certain point, if the bacterial density falls below a critical threshold, the immune system can eliminate the bacteria on its own.

This brings up a thought-provoking question: Why has not our immune system evolved to always have a ε′ value larger than r when activated? While this is complex and cannot be easily answered, excessive inflammation can lead to organ damage or even death. Balancing the immune system is crucial, as emphasized by various studies.

In conclusion, incorporating the host immune system into phage therapy mathematical models allows for a more comprehensive understanding of the complex interactions between phages, bacteria, and the host during treatment. These models can help improve predictions of therapeutic outcomes and guide the development of personalized and effective phage therapy strategies.

OPTIMIZATION OF PHAGE DOSE

Upon a comprehensive scrutiny of diverse elements pertaining to the dynamics among phages, their hosts, and the immune system, our focus now shifts to the determination of the minimum phage dosage essential for successful bacterial eradication, contingent upon immune system functionality. We shall employ an interaction model, tailored from the Leung and Weitz construct [7], to simulate these interactions. This adjusted model accounts for resistant bacterial subpopulations and presumes an innate turnover of baseline immune cells.

(Model XI)

dSdt=r·S·1-BBmax-k·P·S-ε·I·S/(1+BKD)(31) 
dRdt=r·R·1-BBmax-ε·I·R/(1+BKD)(32) 
dPdt=b-1·k·P·S-w·P(33) 
dIdt=α·I·1-IKI·BB+KN+δ·(I0-I)(34) 
B=S+R(35) 

Herein, S and R signify phage-sensitive and resistant bacterial populations, respectively. δ represents the immune cell elimination rate constant, with I0 as the baseline immune cell concentration. We postulate that at the outset, a fraction 0 < φ < 1 of the total bacterial population exhibits resistance. Remaining parameters align with those defined in the sections describing Eq. (27)-(28).

Model simulations reveal a pivotal phage dose threshold beneath which therapeutic endeavors are ineffectual. A diminished efficacy in immune system activation correspondingly elevates this critical dose, as visualized in Fig. 7A in comparison to Fig. 7B. This corroborates previous studies [5, 24] that indicated diminished therapeutic efficacy of phage interventions in neutropenic hosts.

Figure 7
Dose-response simulations based on Model XI with the following default parameter values: B(0) = 106 CFU mL−1, φ = 0.1, r = 1 h−1,k = 10−8 mL CFU−1 h−1,b = 100, w =1 h−1,I0 = 106 cells mL−1,KD = 106 CFU mL−1,KN = 105 CFU mL−1,KI = 107 cells mL−1, ε = 10−6.5 mL CFU−1 h−1,δ = 0.1 h−1, and Bmax = 109 CFU mL−1. (A) and (B) correspond to results assuming α = 1 h−1 and α = 0.85 h−1, respectively.

Furthermore, phage pharmacokinetic dynamics critically influence therapeutic outcomes. Fig. 8 delineates simulation outcomes resulting from modulation of the phage elimination rate constant (w) at doses of 106 PFU mL−1 and 108 PFU mL−1. Remarkably, accelerated phage elimination correlates directly with therapeutic failures.

Figure 8
Effect of phage PK based on Model XI with the following default parameter values: B(0) = 106 CFU mL−1, φ = 0.1, r = 1 h−1,k = 10−8 mL CFU−1 h−1,b = 100, I0 = 106 cells mL−1,KD = 106 CFU mL−1,KN = 105 CFU mL−1,KI = 107 cells mL−1, ε = 10−6.5 mL CFU−1 h−1,δ = 0.1 h−1, and Bmax = 109 CFU mL−1. (A) and (B) correspond to results assuming phage dose of 106 PFU mL−1 and 108 PFU mL−1, respectively.

Central to the success of phage therapy lies the phage’s intrinsic capability to swiftly curtail bacterial concentrations beneath the crucial limit susceptible to immune system regulation. To ascertain rapid bacterial elimination, an adequately high phage dose becomes imperative. Thus, a meticulous exploration of phage pharmacokinetics becomes indispensable in prognosticating therapeutic outcomes.

PHAGE PHARMACOKINETICS

Studies regarding absorption, distribution, metabolism, and excretion (ADME) and dose-ranging have been carried out in animal models; however, a consensus on the optimal route of administration remains elusive and seems to depend heavily on the specific disease being treated. Theoretically, the route should allow phages to proliferate sufficiently to attain the inundation threshold in the target organ.

Research indicates the presence of non-linear pharmacokinetics in phages. While Chow et al. [25] found a dose-proportional surge in phage exposure with doses of 107 and 109 PFU/head administered intratracheally in healthy mice, a separate study conducted by Lin et al. [26] demonstrated diminishing phage clearance with escalating doses when administered intravenously to rats at levels of 106, 109, and 1011 PFU/head. Phages are reportedly eliminated by the innate immune system through the reticuloendothelial system [27], a process that potentially gives rise to nonlinear pharmacokinetics due to the system’s finite capacity.

Predominantly, research posits a rapid clearance of phages from plasma within a 48-hour frame, with a preferential distribution to the liver and spleen, followed by a lesser extent to the lungs. It is plausible that systemic inflammation facilitates accelerated phage clearance [27]. Personal accounts from preclinical developments affirm the non-linear pharmacokinetics in phages, highlighting a considerable variation in pharmacokinetics between different phages. It becomes imperative, therefore, to engage in detailed studies to understand the variable pharmacokinetic behaviors across different phage entities.

THE NECESSITY FOR ADMINISTERING MULTIPLE DOSES

Theoretically speaking, the administration of successive doses in active phage therapy is deemed unnecessary due to the replicative nature of phages, a concept often described with the term “autodosing.” Nonetheless, a considerable number of past clinical case studies have documented the utilization of multiple phage doses [17, 19, 20, 28], a strategy presumably derived from the conventional belief that repeated doses could potentially enhance drug exposure.

However, simulations indicate that in scenarios involving active therapy with phage doses below the inundation threshold, the strategy of administering multiple doses might indeed be detrimental (Fig. 9A). While intermittent phage dosing diminishes host bacterial density, it fails to achieve total eradication. This results in diminished phage multiplication owing to periodic bacterial reduction, an outcome that negatively impacts the ultimate clinical result. Conversely, in passive therapy, employing recurrent administration of doses exceeding the inundation threshold is anticipated to function as presumed, where increased frequency of administration facilitates greater suppression of bacterial populations (Fig. 9C).

Figure 9
Simulations utilizing Model I to investigate the consequences of administering repeat dose at concentrations of 106 PFU mL−1 (panels A and B) and 108 CFU mL−1 (panels C and D) at 24 hours (panels A and C). In the context of an active therapy regimen where doses below the inundation threshold are applied, supplementary dosing results in diminished phage multiplication and subsequently higher nadir host densities. Contrastingly, in a passive therapy scenario where dosages exceed the inundation threshold, an additional dose instantaneously and effectively suppresses host density upon administration. The parameters used for the simulations are: B(0) = 106 CFU mL−1,r = 1 h−1,k = 10−7 mL CFU−1 h−1,b = 100, w = 0.5 h−1, and Bmax = 1010 CFU mL−1. Panels A and B represent outcomes at phage doses of 106 PFU mL−1 and 108 PFU mL−1, respectively.

These findings emphasize the critical role of selecting an appropriate dosing regimen grounded on dose potency. Generally, passive therapy offers a more straightforward treatment strategy, permitting more predictable outcomes without the necessity to factor in the phage proliferation threshold. However, delineating between active and passive therapy is not always straightforward, as it is predicated upon the comparative magnitude of phage concentration against the inundation threshold at the target site, which establishes the therapy’s active or passive nature. Ultimately, the choice to employ multiple dosages should invariably be founded on meticulous examination of the phage’s PKPD characteristics.

ADVANCED MODELING TECHNIQUES AND FUTURE DIRECTIONS IN PHAGE THERAPY RESEARCH

In this last section, we will discuss advanced modeling techniques and future directions in phage therapy research, including the integration of spatial effects, stochasticity, and multi-scale modeling approaches. We will also explore how these advanced techniques can further improve our understanding of phage-bacteria interactions and guide the development of more effective phage therapy strategies.

Spatial effects in phage-bacteria interactions

While the models discussed so far have largely focused on well-mixed systems, in reality, phage-bacteria interactions often occur in spatially structured environments, such as biofilms or host tissues. Spatial heterogeneity can have a significant impact on phage-bacteria dynamics, affecting the spread of phages, the development of bacterial resistance, and the overall success of phage therapy. To incorporate spatial effects, researchers can use techniques such as:

  • 1. Partial differential equations (PDEs): PDE models can capture the spatial distribution of phages and bacteria by adding diffusion terms to the basic model equations, allowing researchers to study the impact of spatial heterogeneity on phage-bacteria interactions. PDE models can be implemented using various software libraries in different programming languages – such as py-pde in Python or MATLAB.

  • 2. Agent-based models (ABMs): ABMs are computational models that represent individual phages, bacteria, and host cells as discrete agents, allowing researchers to simulate spatially explicit interactions between these entities and explore the effects of local interactions on the dynamics of the system. NetLogo (https://ccl.northwestern.edu/netlogo/) is a free software that provides a user-friendly interface for creating ABMs.

Stochastic effects and intrinsic noise

Real-world phage-bacteria systems are often subject to stochastic effects due to the inherent randomness of biological processes, such as phage adsorption, bacterial replication, or the host immune response. To account for stochastic effects, researchers can employ techniques such as:

  • 1. Stochastic differential equations (SDEs): SDEs incorporate noise terms into the basic model equations to capture the random fluctuations in phage and bacterial populations, allowing researchers to study the impact of stochastic effects on phage-bacteria dynamics.

  • 2. Gillespie algorithm: The Gillespie algorithm is a widely used stochastic simulation algorithm that can model the discrete, random nature of molecular interactions in phage-bacteria systems, providing a more realistic representation of the underlying processes.

Multi-scale modeling

Phage-bacteria interactions occur at various scales, from the molecular and cellular levels to the population and community levels. To capture the complexity of these interactions across multiple scales, researchers can develop multi-scale models that integrate processes at different levels of organization. Examples of multi-scale modeling approaches include:

  • 1. Coupling ODE/PDE models with ABMs: Researchers can combine ODE/PDE models, which describe population-level dynamics, with ABMs, which capture individual-level interactions, to create a more comprehensive representation of phage-bacteria systems.

  • 2. Hierarchical models: Hierarchical models can be used to represent the relationships between processes occurring at different scales, such as the impact of phage-bacteria interactions on host tissue dynamics or the role of the microbiome in shaping phage therapy outcomes.

Future directions in phage therapy research

The continued development and application of advanced modeling techniques will undoubtedly lead to new insights into the complex dynamics of phage-bacteria systems and help guide the development of more effective phage therapy strategies. Some promising future directions in phage therapy research include:

  • 1. Personalized phage therapy: Developing models that incorporate patient-specific data, such as the composition of their microbiome or the characteristics of their immune system, to design personalized phage therapy strategies.

  • 2. Synthetic biology and engineered phages: Exploring the use of synthetic biology techniques to engineer phages with improved therapeutic properties, such as enhanced bacterial specificity, reduced immunogenicity, or the ability to overcome bacterial resistance mechanisms.

  • 3. Phage-bacteria coevolution: Studying the coevolutionary dynamics of phages and bacteria, which can influence the long-term success of phage therapy by shaping the emergence of resistance and the evolution of phage infectivity. Developing models that account for these coevolutionary processes can help optimize phage therapy strategies and improve their efficacy over time.

  • 4. Phage cocktails and community-level interactions: Investigating the use of phage cocktails, which are mixtures of multiple phages, to target a broader range of bacterial pathogens and prevent the development of resistance. Developing models that capture the community-level interactions between multiple phages and bacterial species can provide insights into the design of more effective phage cocktails.

  • 5. Integrating experimental and modeling approaches: Combining experimental data with advanced modeling techniques can improve the accuracy and predictive power of phage therapy models. This integration can help bridge the gap between in vitro, in vivo, and clinical studies, ultimately leading to more effective and personalized phage therapy strategies.

In conclusion, advanced modeling techniques and future research directions in phage therapy hold great promise for deepening our understanding of phage-bacteria interactions and improving the effectiveness of phage therapy. By incorporating spatial effects, stochasticity, and multi-scale approaches, researchers can develop more comprehensive and realistic models that better capture the complex dynamics of phage-bacteria systems. These models can guide the development of novel phage therapy strategies, including personalized treatments, engineered phages, and phage cocktails, ultimately contributing to the successful application of phage therapy in clinical settings.

DISCUSSION

The intricate dynamics between phages, host bacteria, and the immune system have long held scientific intrigue. As this review demonstrates, understanding these interactions can significantly influence the application and success rate of phage therapy, especially in the face of bacterial resistance. Mathematical models offer valuable insights into the mechanisms underpinning the interplay among these systems, underscoring the importance of meticulous PKPD evaluations in predicting therapeutic outcomes.

In particular, this article focused on the pivotal role that the immune system plays in determining the efficacy of phage therapy. Model simulations accentuated a critical phage dose threshold, beneath which therapeutic endeavors are rendered ineffectual. This threshold’s sensitivity to immune system efficiency, as depicted in Fig. 7A and B, echoes previous reports which alluded to the diminished therapeutic efficacy in neutropenic hosts. The immune system’s capacity to effectively respond and cooperate with phage interventions has proven paramount in achieving favorable outcomes.

This review also accentuated the importance of the pharmacokinetic characteristics of phages, unveiling that the rapidity of phage elimination directly correlates with the outcome of the therapy. Fig. 8, for instance, elucidates that expediting phage elimination corresponds to therapeutic failures. This revelation underscores the necessity for optimized phage dosage, ensuring not only immediate bacterial eradication but also establishing a bacterial concentration conducive to immune system control.

Furthermore, the article provides an enriched perspective on the existing mathematical models. The incorporation of resistant bacterial subpopulations and assumptions on baseline immune cell turnover in the modified model, adapted from Leung and Weitz [7], enhances the robustness and applicability of the model in real-world scenarios. Such refinements in mathematical constructs propel our comprehension of phage dynamics and can significantly guide future research and therapeutic designs.

In conclusion, phage therapy, despite its promising prospects, is riddled with complexities, most notably its interplay with the host’s immune system and bacterial resistance mechanisms. This review has charted pivotal thresholds and dynamics that inform the application of this therapy, reinforcing the necessity for rigorous PKPD studies and careful therapeutic designs. As bacterial resistance to antibiotics intensifies globally, ensuring the efficacy of alternative treatments like phage therapy becomes even more crucial. Mathematical modeling of phage PKPD not only contribute to the academic discourse but also provide tangible insights for clinicians and pharmacologists aiming to harness the potential of phage therapy effectively. Future endeavors might focus on expanding the current models to consider a broader spectrum of bacterial infections and host responses, solidifying the foundation for phage therapy’s widespread clinical implementation.

Notes

Conflict of Interest:- Author: Nothing to declare

- Reviewers: Nothing to declare

- Editors: Nothing to declare

Reviewer:This article was invited and reviewed by the editors of TCP.

References

    1. Gordillo Altamirano FL, Barr JJ. Phage Therapy in the Postantibiotic Era. Clin Microbiol Rev 2019;32:e00066-18
    1. Strathdee SA, Hatfull GF, Mutalik VK, Schooley RT. Phage therapy: from biological mechanisms to future directions. Cell 2023;186:17–31.
    1. Loc-Carrillo C, Abedon ST. Pros and cons of phage therapy. Bacteriophage 2011;1:111–114.
    1. Cairns BJ, Timms AR, Jansen VA, Connerton IF, Payne RJ. Quantitative models of in vitro bacteriophage-host dynamics and their application to phage therapy. PLoS Pathog 2009;5:e1000253
    1. Roach DR, Leung CY, Henry M, Morello E, Singh D, Di Santo JP, et al. Synergy between the host immune system and bacteriophage is essential for successful phage therapy against an acute respiratory pathogen. Cell Host Microbe 2017;22:38–47.e4.
    1. Delattre R, Seurat J, Haddad F, Nguyen TT, Gaborieau B, Kane R, et al. Combination of in vivo phage therapy data with in silico model highlights key parameters for pneumonia treatment efficacy. Cell Reports 2022;39:110825
    1. Leung CY, Weitz JS. Modeling the synergistic elimination of bacteria by phage and the innate immune system. J Theor Biol 2017;429:241–252.
    1. Rodriguez-Gonzalez RA, Leung CY, Chan BK, Turner PE, Weitz JS. Quantitative models of phage-antibiotic combination therapy. mSystems 2020;5:e00756-19
    1. Smug BJ, Majkowska-Skrobek G, Drulis-Kawa Z. PhREEPred: phage resistance emergence prediction web tool to foresee encapsulated bacterial escape from phage cocktail treatment. J Mol Biol 2022;434:167670
    1. Nilsson AS. Cocktail, a computer program for modelling bacteriophage infection kinetics. Viruses 2022;14:2483.
    1. Weitz JS, Hartman H, Levin SA. Coevolutionary arms races between bacteria and bacteriophage. Proc Natl Acad Sci U S A 2005;102:9535–9540.
    1. Li G, Leung CY, Wardi Y, Debarbieux L, Weitz JS. Optimizing the timing and composition of therapeutic phage cocktails: a control-theoretic approach. Bull Math Biol 2020;82:75.
    1. Bull JJ, Vegge CS, Schmerer M, Chaudhry WN, Levin BR. Phenotypic resistance and the dynamics of bacterial escape from phage control. PLoS One 2014;9:e94690
    1. Wright A, Hawkins CH, Anggård EE, Harper DR. A controlled clinical trial of a therapeutic bacteriophage preparation in chronic otitis due to antibiotic-resistant Pseudomonas aeruginosa; a preliminary report of efficacy. Clin Otolaryngol 2009;34:349–357.
    1. Sarker SA, Berger B, Deng Y, Kieser S, Foata F, Moine D, et al. Oral application of Escherichia coli bacteriophage: safety tests in healthy and diarrheal children from Bangladesh. Environ Microbiol 2017;19:237–250.
    1. Jault P, Leclerc T, Jennes S, Pirnay JP, Que YA, Resch G, et al. Efficacy and tolerability of a cocktail of bacteriophages to treat burn wounds infected by Pseudomonas aeruginosa (PhagoBurn): a randomised, controlled, double-blind phase 1/2 trial. Lancet Infect Dis 2019;19:35–45.
    1. Duplessis C, Biswas B, Hanisch B, Perkins M, Henry M, Quinones J, et al. Refractory Pseudomonas bacteremia in a 2-year-old sterilized by bacteriophage therapy. J Pediatric Infect Dis Soc 2018;7:253–256.
    1. Khawaldeh A, Morales S, Dillon B, Alavidze Z, Ginn AN, Thomas L, et al. Bacteriophage therapy for refractory Pseudomonas aeruginosa urinary tract infection. J Med Microbiol 2011;60:1697–1700.
    1. LaVergne S, Hamilton T, Biswas B, Kumaraswamy M, Schooley RT, Wooten D. Phage therapy for a multidrug-resistant Acinetobacter baumannii craniectomy site infection. Open Forum Infect Dis 2018;5:ofy064
    1. Schooley RT, Biswas B, Gill JJ, Hernandez-Morales A, Lancaster J, Lessor L, et al. Development and use of personalized bacteriophage-based therapeutic cocktails to treat a patient with a disseminated resistant Acinetobacter baumannii infection. Antimicrob Agents Chemother 2017;61:e00954-17
    1. Kortright KE, Chan BK, Koff JL, Turner PE. Phage therapy: a renewed approach to combat antibiotic-resistant bacteria. Cell Host Microbe 2019;25:219–232.
    1. Payne RJ, Phil D, Jansen VA. Phage therapy: the peculiar kinetics of self-replicating pharmaceuticals. Clin Pharmacol Ther 2000;68:225–230.
    1. Abedon ST, Hyman P, Thomas C. Experimental examination of bacteriophage latent-period evolution as a response to bacterial availability. Appl Environ Microbiol 2003;69:7499–7506.
    1. Tiwari BR, Kim S, Rahman M, Kim J. Antibacterial efficacy of lytic Pseudomonas bacteriophage in normal and neutropenic mice models. J Microbiol 2011;49:994–999.
    1. Chow MY, Chang RY, Li M, Wang Y, Lin Y, Morales S, et al. Pharmacokinetics and time-kill study of inhaled antipseudomonal bacteriophage therapy in mice. Antimicrob Agents Chemother 2020;65:e01470-20
    1. Lin YW, Chang RY, Rao GG, Jermain B, Han ML, Zhao JX, et al. Pharmacokinetics/pharmacodynamics of antipseudomonal bacteriophage therapy in rats: a proof-of-concept study. Clin Microbiol Infect 2020;26:1229–1235.
    1. Hodyra-Stefaniak K, Miernikiewicz P, Drapała J, Drab M, Jończyk-Matysiak E, Lecion D, et al. Mammalian host-versus-phage immune response determines phage fate in vivo . Sci Rep 2015;5:14802.
    1. Nir-Paz R, Gelman D, Khouri A, Sisson BM, Fackler J, Alkalay-Oren S, et al. Successful treatment of antibiotic-resistant, poly-microbial bone infection with bacteriophages and antibiotics combination. Clin Infect Dis 2019;69:2015–2018.

Metrics
Share
Figures

1 / 9

ORCID IDs
PERMALINK