Understanding biofilm--phage interactions in cystic fibrosis patients using mathematical frameworks

Blessing O. Emerenini; Doris Hartung; Ricardo N. G. Reyes Grimaldo; Claire Canner; Maya Williams; Ephraim Agyingi; Robert Osgood

doi:10.1515/cmb-2024-0018

Article Open Access

Understanding biofilm--phage interactions in cystic fibrosis patients using mathematical frameworks

Blessing O. Emerenini , Doris Hartung , Ricardo N. G. Reyes Grimaldo , Claire Canner , Maya Williams , Ephraim Agyingi and Robert Osgood

Published/Copyright: January 31, 2025

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Computational and Mathematical Biophysics Volume 13 Issue 1

Abstract

When planktonic bacteria adhere together to a surface, they begin to form biofilms, communities of bacteria. Biofilm formation in a host can be extremely problematic if left untreated, especially since antibiotics can be ineffective in treating the bacteria. Certain lung diseases such as cystic fibrosis can cause the formation of biofilms in the lungs and can be fatal. With antibiotic-resistant bacteria, the use of phage therapy has been introduced as an alternative or an additive to the use of antibiotics to combat biofilm growth. Phage therapy utilizes bacteriophages that attack bacteria, to penetrate and eradicate biofilms. To evaluate the effectiveness of phage therapy against biofilm bacteria, we create a phage-biofilm model with ordinary differential equations and a stochastic model. We considered three possible cases for the stability of disease-free equilibrium; by model simulations and parameter alterations in both models, we investigated the effect of bacteria–phage interactions alongside biofilm bacteria cell detachment from the biofilm phase to the planktonic phase, and how these will affect the efficiency of phage therapy against bacteria. Our results show that, by increasing the phage mortality rate, the biofilm growth can be balanced, which makes it more vulnerable to antibiotics. Thus, phage therapy is an effective aid in biofilm treatment.

Keywords: biofilm; phages; mathematical modeling; stochastic modeling

MSC 2010: 92-XX; 92D30; 60J28; 34D20

1 Introduction

The presence of pathogenic microorganisms in our environment entails enormous problems for humans and livestock. The problem of pathogenic microorganisms is even grievous when they reside in the host [17]. Bacteria is one such pathogenic microorganism, and they prefer to live in communities called biofilms.

Biofilms are aggregations of bacteria on immersed surfaces and interfaces, in which the cells are embedded in a self-produced layer of extracellular polymeric substances (EPS). The EPS gives them protection against mechanical washout and antibiotics. The formation of biofilm is often considered a virulence factor [32]. Given the role that biofilms play in resistance, new treatments are promoted that aim at penetrating the biofilm matrix and attacking the individual cells in the biofilm. Dissolved growth-limiting substrates such as oxygen diffuse into the biofilm and undergo a reaction with bacteria. In many instances, in well-developed biofilms, such growth-limiting substrates might only be able to penetrate the biofilm over relatively thin active outer layers and substantial inactive inner layers may form. Several studies have shown that there are some immaterial substances and microbes that can also penetrate the biofilm matrix, one such microbes is the bacteriophages.

Bacteriophage, also known informally as phage, is a virus that infects and replicates within bacteria and archaea, and it is among the most common disease entities in the biosphere. Bacteriophages are exclusively used as therapeutic agents to treat infections caused by pathogenic bacteria. Application of phage therapy was dated more than a century ago with poor understanding of its potentials [6,36] and so was overshadowed in western medicine until the emergence of bacterial strains that were resistant to antibiotics [1]. The use of phage in treatments has several potential advantages over the use of antibiotics [12,27].

This study focuses on the interaction between biofilm and bacteriophage in phage therapy. There are two main approaches to studying phage–biofilm interactions: the experimental approach and the mathematical modeling approach. Both methods are widely utilized, and each has advantages; the utilization of mathematical models as a proxy to study disease dynamics is an extremely useful tool when it comes to the observation and prevention of infections in humans [10,26,29,37]. Through the use of mathematical modeling, one is able to implement approaches in which infectious outbreaks can be predicted, assessed, and controlled [10,29,38]. To represent these dynamics, a form of model must be chosen with appropriate assumptions. These mathematical models can range from equation-based modeling, such as ordinary, partial, and stochastic differential equations, to agent-based modeling [19,25]. From different studies that utilized mathematical modeling, findings have been made on in-host disease dynamics. Through the use of differential equations in a study conducted by Beke et al., the disease dynamics in a bacteria–phage interaction are examined. The results of this study gathered that disease dynamics can be contingent on environmental factors such as a change in pH or temperature [3]. In a separate study conducted by Bardina et al., a different approach involving a stochastic model is used to analyze these dynamics. Upon these findings, there existed equilibria such that pathogens could be eliminated from the host or could persist depending on the levels of noise within the environment [2]. Similarly, through differential equations and Monte Carlo simulations, Sinha et al. evaluated which mathematical model is adequate for modeling the dynamics between phage and bacteria. In this study, it was found that disease dynamics can differ if there are spatial restrictions introduced to the model and the type of model used to describe such dynamics should reflect this restriction [40]. These mathematical models are not only used for examining in-host dynamics but also used to address various phenomena. Some of these phenomena include spatial phenomena as previously mentioned, evolutionary game theory, and dynamic optimization.

Mathematical models for bacterial biofilm over the years have greatly helped in the understanding of biofilm processes such as biofilm formation and growth, detachment, and its inducers [5,8,15,16,23]. Many of the experimental and modeling studies of biofilm–phage interactions and their interplay have focused on biofilm formed on surfaces other than in-host, mathematical modeling that focuses on the biofilm–phage interactions, and interplay in immunocompromised patients is still in its infancy. Several immunocompromised patients suffer from biofilm infection, and we will be considering the case of biofilm formation in the lungs of cystic fibrosis patients.

Over the last two decades, a large number of models have been produced to represent the interactions between bacteria, biofilm, and bacteriophages. Some of these models involve ordinary differential equations (ODEs), partial differential equations (PDEs), agent-based modeling, or stochastic models, but all of which attempt to recreate results that could be produced in an experiment to better and more quickly understand and predict these interactions. In a recent review from Sinha et al., we find models utilizing PDEs, ODEs, stochastic differential equations, and Monte Carlo simulations. This review compares models subject to spatial constraints with those without spatial constraints [40].

In this study, the objective is to develop a mathematical framework to understand the different factors that contribute most during bacterial–phage interactions in biofilm setting and planktonic phase.

2 Methods

2.1 Basic model assumptions

We develop a mathematical model that describes the dynamics of biofilm–phage interactions in matured biofilms, under the condition that both bacterial and phage populations are high enough; the interaction between bacteria, phages, and nutrients does not depend on their location; the biofilm-planktonic system is closed; and nutrients enter the system at a predetermined rate ( f ( n ) ). There are two specific regions involved in this system, namely, the biofilm phase, which is the region where the bacterial cells accumulate and form biofilm; and the Planktonic phase, which is the flow region of the bronchiole comprising of air and fluid (Figure 1). We assume that (a) there are two different forms of bacteria–phage interaction taking place in the system: one interaction leads to infection and ultimately to the replication of phages; the other interaction leads to the conversion of bacteria or phages from one region to another [16,18], (b) conversion of biofilm bacteria to planktonic bacteria is induced by phages, (c) phages conversion rate between biofilm and the planktonic phase is assumed to be constant ( p and q ), and (d) the conversion of phages from one region to the other does not change their characteristics.

Figure 1

Schematic of biofilm in bronchiole. This is a schematic representation of the biofilm that formed on the wall of the bronchiole, showing the biofilm region and the planktonic region. Source: This figure has been created by the authors.

2.2 Deterministic model – ODE

The model describing the biofilm–phage interactions is formulated as a deterministic model of a system of six ordinary differential equations. The dependent variables are n , B , P , V B , V P and I . At a time t , the variable n denotes the concentration of the bacteria growth limiting nutrient substrate, B denotes the bacteria cells in the biofilm region while P denotes the bacterial cells in the planktonic region, V B and V P denotes the viral load of bacteriophages in the biofilm and the planktonic regions respectively; I is the concentration of all the infected bacteria cells from the planktonic and biofilm regions. The dynamics for substrate nutrients are established by the functional f ( n ) which under certain conditions, defined in the next section, can facilitate desired scenarios like nutrient saturation. The bacterial population growth is based on a monod function of n and decreases from their intrinsic mortality rates, conversion between the biofilm and planktonic regions, or infection as a consequence of the interaction with the corresponding phage population. The model captures the detachment of bacteria cells from the biofilm induced by bacteriophages and reattachment of bacteria to the biofilm, and this meshes well with the life cycle of biofilms.

The governing equations for our deterministic model are given by the following system:

(1) d d t n = f ( n ) − ( λ B B + λ P P ) n n + k

(2) d d t B = λ B B n n + k − ϕ 1 B V B − γ 1 V B ζ 1 + V B B + γ 2 B ζ 2 + B P − μ B B

(3) d d t P = λ P P n n + k − ϕ 2 P V P + γ 1 V B ζ 1 + V B B − γ 2 B ζ 2 + B P − μ P P

(4) d d t V B = β ϕ 1 V B B − c 1 V B − q V B + p V P

(5) d d t V P = β ϕ 2 V P P − c 2 V P + q V B − p V P

(6) d d t I = ϕ 1 V B B + ϕ 2 V P P − 10 τ n n + k I − α I .

Equations (2) and (3) describe the bacterial growth in the planktonic and biofilm regions, the equations also describe predation of bacteria by the phages and biofilm cell detachment which forms a direct coupling of the system. Equations (4) and (5) describe the phage growth in the biofilm and planktonic region, while equation (6) keeps track of all infected bacterial cells in both the biofilm and planktonic regions. Equation (1) describes the consumption of nutrient by B and P . The flow diagram of the model is presented in Figure 2. The parameters of this model are presented in Table 1, and one of the parameters of interest is the burst size that determines the average number of phage releases per bacterium and could vary from 10 to 100 for DNA transducing bacteriophages to about 20,000 pfu for the RNA viruses.

Figure 2

Flow diagram of the deterministic ODE model. This is a schematic representation of the deterministic ODE model, showing the flows and connections with the parameters. These parameters are also presented in Table 1 with the actual descriptions of the parameter, the values, and sources. Because of the lack of a better word, we have used “migration” to capture the conversion from biofilm to planktonic region, and this does not imply any spatial component. Source: This figure has been created by the authors.

Table 1

Table of parameters

Symbol	Description	Value	Units	Source
n 0	Initial nutrient concentration	Variable	g m ‒ 3	[16,18]
B 0	Initial biofilm bacteria	Variable	g m ‒ 3	[16,18]
P 0	Initial planktonic bacteria	Variable	g m ‒ 3	[18]
V b 0	Initial biofilm phages	Variable	g m ‒ 3	Assumed
V p 0	Initial planktonic phages	Variable	g m ‒ 3	[18]
I 0	Initial Infected cells	Variable	g m ‒ 3	[18]
p	Phage detachement rate	0.1	d ‒ 1	Assumed
q	Phage reattach rate	0.5	d ‒ 1	Assumed
λ B	biofilm bacteria growth rate	6.0	d ‒ 1	[18]
λ P	Planktonic bacteria growth rate	6.0	d ‒ 1	[18]
τ	Average latency time	0.5	h	[11,18]
k	Monod constant	4.0	g m ‒ 3	[18]
α	Infection decay rate	0.2	d ‒ 1	Assumed
β	Burst size	100	—	[18]
γ 1	Phage induced detach rate	0.6	d ‒ 1	[16]
γ 2	Natural detach rate	0.3	d ‒ 1	[16]
ϕ 1	Predation rate in biofilm	1 0 ‒ 8	m 3 g ‒ 1 d ‒ 1	Assumed
ϕ 2	Predation rate in planktonic	1 0 ‒ 6	m 3 g ‒ 1 d ‒ 1	Assumed
ζ 1	Monod saturation	1 0 2	g m ‒ 3	Assumed
ζ 2	Monod saturation	1 0 4	g m ‒ 3	Assumed
c 1	Phage mortality rate in biofilm	2.1	d ‒ 1	Assumed
c 2	Phage mortality rate in planktonic	2.1	d ‒ 1	Assumed
μ B	Bacteria mortality rate in Biofilm	0.1	d ‒ 1	Assumed
μ P	Bacteria mortality rate in planktonic	0.1	d ‒ 1	Assumed

2.2.1 Basic reproduction number

The system (1)–(6) represents a nonlinear system of ODEs with the interaction of phages and bacteria. We aim to determine the stability of equilibria points in the model through the basic reproduction number, which will describe the expected number of infected bacteria generated by one case of phage infection, using the next generation matrix method as described in previous studies [41,42].

Remark 1

(Assumption for the nutrient rate) The function f ( n ) , for the nutrient rate of change, is such that

f is continuous and differentiable on the interval [ 0 , ∞ ) .
f ( 0 ) = 0 (thus making the origin a steady state) and f ′ ( 0 ) > 0 .
For some fixed n , say x ∗ > 0 , the conditions f ( x ∗ ) = 0 and f ′ ( x ∗ ) < 0 are held.

There are several functions that satisfy the aforementioned conditions such as f ( n ) = n ( 1 − n ) .

Considering the aforementioned properties and system (1)–(6), we compute the Jacobian matrix, which is given by:

(7) J ( X , P ) = f ′ ( n ) − k ( λ B B + λ P P ) ( n + k ) 2 − n λ B n + k − n λ P n + k 0 0 0 k B λ B ( n + k ) 2 † B γ 2 ζ 2 + B − ϕ 1 B − ζ 1 γ 1 B ( ζ 1 + V B ) 2 0 0 λ P P k ( n + k ) 2 γ 1 V B ζ 1 + V B − γ 2 ζ 2 P ( ζ 2 + B ) 2 ‡ γ 1 ζ 1 B ( ζ 1 + V B ) 2 − ϕ 2 P 0 0 β ϕ 1 V B 0 β ϕ 1 B − c 1 − q p 0 0 0 β ϕ 2 V P q β ϕ 2 P − p − c 2 0 − 10 k I τ ( n + k ) 2 ϕ 1 V B ϕ 2 V P ϕ 1 B ϕ 2 P − 10 n τ ( n + k ) − α ,

where

† = n λ B n + k − ϕ 1 V B − γ 1 V B ζ 1 + V B + γ 2 ζ 2 P ( ζ 2 + B ) 2 − μ B ‡ = λ P n n + k − ϕ 2 V P − γ 2 B ζ 2 + B − μ P .

The equilibria points are determined by the zeros of the systems (1)–(6). Considering the disease-free equilibria, which are mainly determined in the absence of a pathogen scenario, that is, when neither phages nor infected bacterial cells are present. From this context, we let I ∗ = 0 , V P ∗ = 0 , and V B ∗ = 0 ; notice this condition causes equations (4)–(6) to be equal to zero, thus reducing the equilibria problem to just find the zeros of the following system

(8) d d t n = f ( n ) − ( λ B B + λ P P ) n n + k d d t B = B λ B n n + k + γ 2 P ζ 2 + B − μ B d d t P = P λ P n n + k − γ 2 B ζ 2 + B − μ P .

Through a quick inspection, we have that the equilibria points of (8) that are biologically relevant are defined by

DFE 1 = ( n ∗ , 0 , 0 , 0 , 0 , 0 ) ; where n ∗ is a zero of the function f DFE 2 = k μ B λ B − μ B , f k μ B λ B − μ B k μ B λ B − μ B + k λ B k μ B λ B − μ B , 0 , 0 , 0 , 0 = n B ∗ , f ( n B ∗ ) μ B , 0 , 0 , 0 , 0 DFE 3 = k μ P λ P − μ P , 0 , f k μ P λ P − μ P k μ P λ P − μ P + k λ P k μ P λ P − μ P , 0 , 0 , 0 = n P ∗ , 0 , f ( n P ∗ ) μ P , 0 , 0 , 0 DFE 4 = ( n ∗ , B ∗ , P ∗ , 0 , 0 , 0 ) described in the following section ,

where n B ∗ and n P ∗ need to be positive, and such that when f is evaluated at those values, f ( n P ∗ ) and f ( n B ∗ ) are also positive, to guarantee biologically relevant scenarios.

We seek to address the stability of each of the equilibria above and under which conditions are each locally asymptotically stable. First, by evaluating the Jacobian (7) at DFE 1 we have that

(9) J 0 , 1 = f ′ ( n ∗ ) − λ B n ∗ n ∗ + k − λ P n ∗ n ∗ + k 0 0 0 0 λ B n ∗ n ∗ + k − μ B 0 0 0 0 0 0 λ P n ∗ n ∗ + k − μ P 0 0 0 0 0 0 − c 1 − q p 0 0 0 0 q − c 2 − p 0 0 0 0 0 0 − 10 n ∗ τ ( n ∗ + k ) − α ,

thus allowing us to obtain the corresponding basic reproduction number for this equilibrium point through the next generation method [42], defining:

(10) R 0 , 1 = max λ B n ∗ μ B ( n ∗ + k ) , λ P n ∗ μ P ( n ∗ + k ) , ± p q ( c 1 + q ) ( c 2 + p ) .

Theorem 1

The disease-free equilibrium 1 ( DFE 1 ) is locally asymptotically stable if n ≠ 0 and R 0 , 1 < 1 .

Proof

Suppose n ∗ = 0 is a zero of f , considering Remark 1 (assumption for the nutrient rate), and in this case, the system naturally becomes unstable since (9) has positive eigenvalues, so it suffices that n * > 0 be a zero of f .

Remains to show R 0 , 1 < 1 for stability: For this, the overall stability of DFE 1 = ( n * , 0 , 0 , 0 , 0 , 0 ) is given by the roots of the characteristic polynomial of (9), which is

det ( s I − J 0 , 1 ) = det s − f ′ ( n ∗ ) λ B n ∗ n ∗ + k λ P n ∗ n ∗ + k 0 0 0 0 s + μ B − λ B n ∗ n ∗ + k 0 0 0 0 0 0 s − λ P n ∗ n ∗ + k + μ P 0 0 0 0 0 0 s + c 1 + q − p 0 0 0 0 − q s + c 2 + p 0 0 0 0 0 0 s + 10 n ∗ τ ( n ∗ + k ) + α = ( s − f ′ ( n * ) ) s + μ B − λ B n * n * + k s − λ P n * n * + k + μ P s + 10 n * τ ( n * + k ) + α × ( s 2 + ( c 1 + q + c 2 + p ) s + ( c 1 + q ) ( c 2 + p ) − p q ) ,

whose roots are all negative if the following holds:

μ B > λ B n * n * + k , λ P n * n * + k < μ P , [ c 1 + q ] [ c 2 + p ] > p q .

Notice that when R 0 , 1 < 1 , all the conditions above hold. Hence, DFE 1 is asymptotically stable. Furthermore, notice that the last inequality arises from the constant term in the quadratic at the end of the characteristic polynomial of (9), as it makes it positive thus having a quadratic polynomial with no change signs, by Decarte’s sign rule such polynomial has no positive real roots, and when the roots are complex they have a negative real part.□

Second for DFE 2 , its corresponding Jacobian matrix is defined by

J 0 , 2 = f ′ ( n B ∗ ) − k ⋆ − q ⋆ λ B − q ⋆ λ P 0 0 0 k ⋆ q ⋆ λ B − μ B s ⋆ − t ⋆ ϕ 1 − γ 1 t ⋆ ζ 1 0 0 0 0 q ⋆ λ P − s ⋆ μ P γ 1 t ⋆ ζ 1 0 0 0 0 0 β t ⋆ ϕ 1 − c 1 − q p 0 0 0 0 q − p − c 2 0 0 0 0 t ⋆ ϕ 1 0 − 10 q ⋆ τ − α ,

where

k ⋆ = k λ B f ( n B ∗ ) μ B ( n B ∗ + k ) 2 , q ⋆ = n B ∗ n B ∗ + k , s ⋆ = f ( n B ∗ ) γ 2 μ B ζ 2 + f ( n B ∗ ) , t ⋆ = f ( n B ∗ ) μ B .

We compute the corresponding basic reproduction number using the next generation matrix method [42], thus:

(11) R 0 , 2 = max k λ B f ( n B * ) ∣ f ′ ( n B * ) ∣ μ B ( n B * + k ) 2 + k λ B f ( n B * ) , ( μ B ξ 2 + f ( n B * ) ) λ P η B * ( n B * + k ) { γ 2 f ( n B * ) + μ P ( μ B ξ 2 + f ( n B * ) ) } , ϕ 1 f ( n B * ) μ B ( c 1 + q ) .

Theorem 2

The disease-free equilibrium 2 ( DFE 2 ) is locally asymptotically stable if n B ∗ ≠ 0 , f ′ ( n B ∗ ) ≫ 0 and R 0 , 2 < 1 .

Proof

Suppose n B * = 0 is a zero of f , considering Remark 1 (assumption for the nutrient rate), and in this case, the system naturally becomes unstable since the DFE 2 Jacobian matrix has positive eigenvalues, so it suffices that n * > 0 be a zero of f .

Remains to show R 0 , 2 < 1 for stability: For this, the overall stability of DFE 2 = n B ∗ , f ( n B ∗ ) μ B , 0 , 0 , 0 , 0 is given by the roots of the characteristic polynomial of DFE 2 Jacobian matrix, which is

det ( s I − J 0 , 2 ) = det s − f ′ ( n B ∗ ) + k ⋆ q ⋆ λ B q ⋆ λ P 0 0 0 − k ⋆ s − q ⋆ λ B + μ B − s ⋆ t ⋆ ϕ 1 + γ 1 t ⋆ ζ 1 0 0 0 0 s − q ⋆ λ P + s ⋆ μ P γ 1 t ⋆ ζ 1 0 0 0 0 0 s − β t ⋆ ϕ 1 + c 1 + q p 0 0 0 0 q s + p + c 2 0 0 0 0 − t ⋆ ϕ 1 0 s + 10 q ⋆ τ + α ,

where

k ⋆ = k λ B f ( n B ∗ ) μ B ( n B ∗ + k ) 2 , q ⋆ = n B ∗ n B ∗ + k , s ⋆ = f ( n B ∗ ) γ 2 μ B ζ 2 + f ( n B ∗ ) , t ⋆ = f ( n B ∗ ) μ B

whose roots are all negative if the following holds:

If f ′ ( n B ∗ ) + n B ∗ λ B n B ∗ + k < k λ B f ( n B ∗ ) μ B ( n B ∗ + k ) 2 + μ B ,
- – Since phages are gone, and bacteria are growing, we expect the mortality rate of bacteria to be negligible when compared with the growth rate.
If n B ∗ λ B n B ∗ + k + k λ B f ( n B ∗ ) f ′ ( n B ∗ ) ( n B ∗ + k ) 2 > μ B
- – This means that the bacteria growth in the biofilm is greater than the bacteria death, and this sounds right since the DFE of phages in the biofilm implies losing more phages and having more bacteria.
If ( c 2 + p ) c 1 + q c 2 > β ϕ 1 f ( n B ∗ ) ( c 2 + p ) μ P
- – This means that the phage loss in the biofilm is greater than phage growth.

Notice that when R 0 , 2 < 1 , the aforementioned conditions hold. Hence, DFE 2 is locally asymptotically stable when R 0 , 2 . This concludes the proof.□

Finally, for DFE 3 , the corresponding Jacobian matrix is given by

J 0 , 3 = f ′ ( n P ∗ ) − k λ P f ( n P ∗ ) μ P ( n P ∗ + k ) 2 − n P ∗ λ B n P ∗ + k − n P ∗ λ P n P ∗ + k 0 0 0 0 n P ∗ λ B n P ∗ + k + γ 2 f ( n P ∗ ) μ P ζ 2 − μ B 0 0 0 0 λ P k f ( n P ∗ ) μ P ( n P ∗ + k ) 2 − γ 2 f ( n P ∗ ) μ P ζ 2 λ P n P ∗ n P ∗ + k − μ P 0 − ϕ 2 f ( n P ∗ ) μ P 0 0 0 0 − c 1 − q p 0 0 0 0 q β ϕ 2 f ( n P ∗ ) μ P − p − c 2 0 0 0 0 0 ϕ 2 f ( n P ∗ ) μ P − 10 n P ∗ τ ( n P ∗ + k ) − α ,

we compute the basic reproduction number using the next generation method [42], thus:

(12) R 0 , 3 = max n P ∗ λ B μ B ( n P ∗ + k ) + γ 2 f ( n P ∗ ) μ P μ B ζ 2 , λ P n P ∗ μ P ( n P ∗ + k ) 1 − k λ P f ( n P ∗ ) ∣ f ′ ( n P ∗ ) ∣ μ P ( n P ∗ + k ) 2 + k λ P f ( n P ∗ ) , β ϕ 2 f ( n P ∗ ) 2 μ P ( c 2 + p ) ± β ϕ 2 f ( n P ∗ ) 2 μ P ( c 2 + p ) 2 + p q ( c 1 + q ) ( c 2 + p ) .

Theorem 3

The disease-free equilibrium 3 ( DFE 3 ) is locally asymptotically stable if n p * ≠ 0 , f ′ ( n p * ) ≫ 0 and R 0 , 3 < 1 .

Proof

Suppose n p * = 0 is a zero of f , considering Remark 1 (assumption for the nutrient rate), in this case, the system naturally becomes unstable since the DFE 3 Jacobian matrix has positive eigenvalues, so it suffices that n p * > 0 be a zero of f .

Remains to show R 0 , 3 < 1 for stability: For this, the overall stability of DFE 3 = n P ∗ , 0 , f ( n P ∗ ) μ P , 0 , 0 , 0 is given by the roots of the characteristic polynomial of the Jacobian matrix J 0 , 3 . The roots are indeed all negative if the following inequalities are satisfied, which have equivalent meaning as those described for DFE 2 :

n P ∗ λ B n P ∗ + k + γ 2 f ( n P ∗ ) μ P ζ 2 < μ B
f ′ ( n P ∗ ) + λ P n P ∗ n P ∗ + k < μ P + k λ P f ( n P ∗ ) μ P ( n P ∗ + k ) 2 and f ′ ( n P ∗ ) λ P n P ∗ n P ∗ + k + k λ P f ( n P ∗ ) ( n P ∗ + k ) 2 > μ P f ′ ( n P ∗ )
β ϕ 2 f ( n P ∗ ) μ P < p + c 2 + q + c 1 and β ϕ 2 f ( n P ∗ ) μ P < p c 1 c 1 + q + c 2 .

When R 0 , 3 < 1 , all these conditions hold. Hence, DFE 3 is locally asymptotically stable when R 0 , 3 < 1 , concluding this proof.□

2.3 Co-existence disease-free equilibrium

The remaining piece is DFE 4 , which is a disease-free equilibrium where there is coexistence of the bacterial populations within the biofilm and the planktonic region. Hence, we need to solve the following system of equations:

(13) 0 = f ( n ) − ( λ B B + λ P P ) n n + k

(14) 0 = λ B n n + k + γ 2 P ζ 2 + B − μ B

(15) 0 = λ P n n + k − γ 2 B ζ 2 + B − μ P .

From (15), we find that the second entry of this equilibrium point holds

(16) B ∗ = ζ 2 λ P n n + k − μ P γ 2 + μ P − λ P n n + k .

From (14), we find that the third entry of this equilibrium point holds

(17) P ∗ = ζ 2 + B ∗ γ 2 μ B − λ B n n + k .

The first entry of this equilibrium point is the solution of the following nonlinear equation:

f ( n ∗ ) = ( λ B B ∗ + λ P P ∗ ) n ∗ n ∗ + k ,

by substituting (16) and (17) in equation results in

(18) f ( n ∗ ) = B ∗ λ B + λ P γ 2 μ B ( n ∗ + k ) − λ B n ∗ n ∗ + k + λ P ζ 2 γ 2 μ B ( n ∗ + k ) − λ B n ∗ n + k n ∗ n ∗ + k .

2.3.1 Biological relevance of coexistence equilibrium

One important characteristic needed of this equilibrium point is its biological relevance, that is, for all three entries to be positive. Notice that this holds if the following scenarios hold

B ∗ > 0 if and only if
- λ P n ∗ > μ P > 0 and ( γ 2 + μ P ) ( n ∗ + k ) > λ P n ∗ > 0
- OR
- μ P > λ P n ∗ > 0 and λ P n ∗ > ( γ 2 + μ P ) ( n ∗ + k ) > 0
P ∗ > 0 if and only if
- B ∗ > 0 and μ B > λ B n ∗ n ∗ + k .

Furthermore, notice that the second scenario for the positivity of B ∗ mentioned earlier is unlikely to occur since it requires that a linear re-scaling of μ P would be larger than λ P n ∗ with the additional condition of λ P n ∗ > μ P . Moreover, the condition on the positivity of P ∗ is such that the right-hand side of (18) is strictly positive for n ∗ > 0 , and hence, the possible intersection with f ( n ∗ ) is not guaranteed to occur in experimental setting given the following circumstances:

f ( 0 ) = 0 .
For some x ∗ > 0 , the conditions f ( x ∗ ) = 0 and f ′ ( x ∗ ) < 0 are held.
g ( n ∗ ) ≔ B ∗ λ B + λ P γ 2 μ B ( n ∗ + k ) − λ B n ∗ n ∗ + k + λ P ζ 2 γ 2 μ B ( n ∗ + k ) − λ B n ∗ n ∗ + k is such that g ( 0 ) = B ∗ λ B + λ P μ B γ 2 + λ P ζ 2 μ B γ 2 > 0 , has a vertical asymptote at n ∗ = − k , and horizontal asymptote at B ∗ λ B + λ P ( μ B − λ B ) γ 2 + λ P ζ 2 ( μ B − λ B ) γ 2 and a root at n ∗ = − 1 − ( 1 + k ) ( B ∗ + ζ 2 ) μ P μ B B ∗ λ B γ 2 .
g ′ ( n ∗ ) = − B ∗ λ P λ B k γ 2 ( n ∗ + k ) 2 + ζ 2 λ P λ B k γ 2 ( n ∗ + k ) 2 < 0 for all n ∗ , and hence, this is a strictly decreasing function.

In the cases where the maximum values for f ( n ∗ ) are under the asymptote of g ( n ) , we cannot have the conditions for a disease-free equilibrium with bacterial coexistence in an experimental setting and therefore biologically irrelevant.

Remark 2

Notice that among the conditions that we are requiring for B ∗ and P ∗ , we require the absorption of nutrients for both bacteria to be smaller than the respective mortalities, and at the same time, the rate of nutrient absorption is larger than the mortalities in one region but not the other.

3 Numerical simulation

All the computations were done using Matlab, the parameter values of the model equations are listed in Table 1, which shows the parameter values, units, and sources. For the nutrient input function f ( n ) , we fixed a logistic growth given by f ( n ) = n ( 1 − n ⁄ k ) , where k is the Monod half-saturation constant. Due to the novelty of this study, some of the parameter values are assumed based on similar studies such as [16]. In this section, we solve and simulate the deterministic model described earlier to understand the system’s behavior for a given set of known initial conditions; this is compared with the sample paths from the stochastic model, which captures some levels of randomness present in the model through a continuous time Markov chain (CTMC) model.

Parameter sensitivity was performed in the subsequent section to determine which of the parameter values have a strong impact on the model generally. We let the program run for at least 15 units of time to see the model behavior in its completeness.

3.1 ODE model: Phage burst size controls biofilm growth

The first simulation experiments investigates the effect of phage burst size on the biofilm growth. Here, we have considered a situation of a steady supply of nutrients to the biofilm, and we have varied the burst size as 10, 50, 150, and 250. These results are presented in Figure 3, in which, we have considered only biofilm cells and phages in the biofilm and in the planktonic phase, with the same initial conditions in each of the simulations. The biofilm is assumed to be growing already with a moderate viral load of phages at the onset of the simulation. As the simulation progresses, we observed that (i) the phages in the biofilm increases as the biofilm increases, (ii) the increased phage growth in the biofilm did not clear the bacterial cells in the biofilm but rather stabilizes the bacterial growth in the biofilm, and (iii) the planktonic phages grows and balances the bacterial cell growth in the planktonic phase. The rationale here is to ascertain the effect of the biofilm in the case of targeted treatment of injecting the phages into the biofilm directly; so with this, we can compare our results with the experimental data. Several studies show that Planktonic bacterial cells grow more rapidly than bacteria cells within a biofilm; therefore, the phage burst size in the biofilm could be several-fold smaller and the infection cycle takes even longer [7,13,20–22,24,35,39]. Remarkably there is a standard procedure for the determination of phage burst size, which is defined as the number of phage progeny produced per infected bacterial cell [9,11,14,30]. Phage burst size differs from phage to phage depending on the lysis time.

$Figure 3 Variability of phage burst size. Left - The average value of the biofilm cells B B , planktonic cells P P , biofilm phage V B {V}_{B} and the planktonic phage V P {V}_{P} , this is plotted for different values of the burst size β \beta . Right - The average possible value of the biofilm phage for different burst sizes (black solid lines) is compared with the data from [43] for different lysis times for comparable burst sizes (red). Source: This figure has been created by the authors.$

Figure 3

Variability of phage burst size. Left - The average value of the biofilm cells B , planktonic cells P , biofilm phage V B and the planktonic phage V P , this is plotted for different values of the burst size β . Right - The average possible value of the biofilm phage for different burst sizes (black solid lines) is compared with the data from [43] for different lysis times for comparable burst sizes (red). Source: This figure has been created by the authors.

3.2 Stochastic model - CTMC

If the bacteria (or phage) population is sufficiently small, an ordinary differential equation model is not appropriate, and hence, we utilize a CTMC model, which is continuous in time and discrete in the state space in order to study the variability at the initiation of bacteria clearance during phage treatment therapy and, peak level of phage infection (in the phage–bacteria interaction, phages are seen as the pathogen and the bacteria are the susceptible). To make it simple, we use the same notation for the state variables as in the ordinary differential equation. The state variables are discrete random variables, n , B , I , P ∈ { 0 , 1 , 2 , 3 } and t ∈ [ 0 , ∞ ]

To formulate the CTMC, it is necessary to define the infinitesimal transition probabilities that corresponds to each event in the state variables, and this is outlined in Table 2, which consists of 17 distinct events.

Table 2

Transitions and corresponding probabilities in stochastic model

Table of events
Event	Event description	Transitions	Change ( Δ n , Δ B , Δ P , Δ V b , Δ V p , Δ I )	Probability
1	Availability of nutrient	n → n + 1	( 1 , 0 , 0 , 0 , 0 , 0 )	f ( n ) Δ t + o ( Δ t )
2	Nutrient consumption and bacterial growth	n → n ‒ 1	( ‒ 1 , 1 , 0 , 0 , 0 , 0 )	( λ B B ) n n + k Δ t + o ( Δ t )
		B → B + 1
3	Nutrient consumption and bacterial growth	n → n ‒ 1	( ‒ 1 , 0 , 1 , 0 , 0 , 0 )	( λ P P ) n n + k Δ t + o ( Δ t )
		P → P + 1
4	Bacteria migration	B → B ‒ 1	( 0 , ‒ 1 , + 1 , 0 , 0 , 0 )	γ 1 V B ζ 1 + V B B Δ t + o ( Δ t )
		P → P + 1
5	Bacteria migration	P → P ‒ 1	( 0 , 1 , ‒ 1 , 0 , 0 , 0 )	γ 2 B ζ 2 + B P Δ t + o ( Δ t )
		B → B + 1
6	Biofilm bacteria infection by phage	B → B ‒ 1	( 0 , ‒ 1 , 0 , 0 , 0 , 1 )	( ϕ 1 B V B ) Δ t + o ( Δ t )
		I → I + 1
7	Planktonic bacteria infection by phage	P → P ‒ 1	( 0 , 0 , ‒ 1 , 0 , 0 , 1 )	( ϕ 2 P V P ) Δ t + o ( Δ t )
		I → I + 1
8	Biofilm–phage migration	V B → V B ‒ 1	( 0 , 0 , 0 , ‒ 1 , 1 , 0 )	( q V B ) Δ t + o ( Δ t )
		V P → V P + 1
9	Biofilm–phage migration	V P → V P ‒ 1	( 0 , 0 , 0 , 1 , ‒ 1 , 0 )	( p V P ) Δ t + o ( Δ t )
		V B → V B + 1
10	Biofilm–phage gain from infected cells	V B → V B + 1	( 0 , 0 , 0 , 1 , 0 , 0 )	( β ϕ 1 V B B ) Δ t + o ( Δ t )
11	Planktonic–phage gain from infected cells	V P → V P + 1	( 0 , 0 , 0 , 0 , 1 , 0 )	( β ϕ 2 V P P ) Δ t + o ( Δ t )
12	Death of biofilm bacteria	B → B ‒ 1	( 0 , ‒ 1 , 0 , 0 , 0 , 0 )	( μ B B ) Δ t + o ( Δ t )
13	Death of planktonic bacteria	P → P ‒ 1	( 0 , 0 , ‒ 1 , 0 , 0 , 0 )	( μ P P ) Δ t + o ( Δ t )
14	Biofilm–phage death	V B → V B ‒ 1	( 0 , 0 , 0 , ‒ 1 , 0 , 0 )	( c 1 V B B ) Δ t + o ( Δ t )
15	Planktonic–phage death	V P → V P ‒ 1	( 0 , 0 , 0 , 0 , ‒ 1 , 0 )	( c 2 V P P ) Δ t + o ( Δ t )
16	Decay of infected cells	I → I ‒ 1	( 0 , 0 , 0 , 0 , 0 , ‒ 1 )	( 10 n τ ( n + k ) ) Δ t + o ( Δ t )
17	Death of infected cells	I → I ‒ 1	( 0 , 0 , 0 , 0 , 0 , ‒ 1 )	( α I ) Δ t + o ( Δ t )

3.2.1 CTMC analysis

For the CTMCs, we numerically simulate the sample paths to visually observe the stochasticity in the model, we also evaluated the frequency of infection to determine the peak number of infected bacteria and peak phage–bacteria infection. For the sample paths, we simply compare our results with that of the ordinary differential equations.

3.2.2 Sample paths

An example of the sample paths that result from the continuous time Markov chain model is shown in Figure 4, where the sample paths are captured by the blue, orange, purple, and yellow lines, whereas the ODE model is captured by the black line presented in the subplots. We observe that these sample paths generally aligned with the population average response that is captured by the ODE model. The sample paths of the CTMC model show the potential variability in timing of the peak level of infection and the peak number of infected bacteria, which is captured by the infection frequency in Figure 5. Due to limitation in computational memory, the dependent variables are presented as volume fractions with reduced initial values as shown in Table 2.

$Figure 4 Sample paths and deterministic plots. This is the sample paths to the CMTC model (left) and the deterministic model plots (right) with black lines, and due to insufficient memory space, we used reduced initial conditions for all the dependent variables, and volume fractions for the phages. Source: This figure has been created by the authors.$

Figure 4

Sample paths and deterministic plots. This is the sample paths to the CMTC model (left) and the deterministic model plots (right) with black lines, and due to insufficient memory space, we used reduced initial conditions for all the dependent variables, and volume fractions for the phages. Source: This figure has been created by the authors.

Figure 5

Peak infection. This shows the frequency of infection of bacteria up to peak infection (left) and the corresponding time to reach the peak infection (right), which is also captured by the mean time to reach the highest frequency of infection of bacteria. Source: This figure has been created by the authors.

3.2.3 Time to peak infection and peak number of infected biofilm bacteria

We asked whether bacteria infection will reach peak infection in a shorter time; to investigate this, we calculated the mean ( ± SD ) of time to peak infection for the bacteria in the biofilm. This is presented in Figure 5 showing that we can attend peak infection with just one phage in the system within a limited time. By introducing few bacteriophages, we observed a large amount of infected bacterial cells resulting from the interaction, and this shows that there were a large replication of the viruses. Interestingly, this happened within a short period of time. Even though we do not know how long it might take phage therapy to work, experimental data have shown that treatment of bacteria infection could be achieved in a period as short as 10 days and up to 8 weeks [33], and this is in agreement with our finding.

4 Parameter sensitivity analysis

We conduct a sensitivity analysis on the parameter ranges outlined in Table 3 for the ODE models, utilizing a uniform distribution for the parameter values. To achieve this, we apply Latin hypercube sampling (LHS), a technique first introduced by McKay et al. [30,31], in conjunction with the partial rank correlation coefficient (PRCC) as the statistical sensitivity measure. This approach enables a comprehensive exploration of the parameter space defined by the intervals specified in Table 3. Unlike traditional methods that assess one parameter at a time while keeping others fixed at baseline values, the LHS/PRCC method allows for a global exploration of the multidimensional parameter space. LHS is a stratified Monte Carlo sampling technique without replacement, which provides an unbiased estimate of the average model output using a limited number of samples. The PRCC is particularly effective for parameters that exhibit nonlinear and monotonic relationships with the output measure.

Table 3

Baseline parameter values are used in all simulations

	Parameter	Values
		Baseline	Range
λ B	Biofilm bacteria growth rate	0.06
λ P	Planktonic bacteria growth rate	0.6931	(0.4, 0.8)
k	Monod constant	6.3	(2–8)
γ 1	Phage induced detachment rate	0.4	(0.1, 2)
γ 2	Natural detachment rate	0.08	(0.02, 0.14)
ϕ 1	Biofilm predation rate	0.075	(0.05, 0.1)
ϕ 2	Planktonic predation rate	0.075	(0.05, 0.1)
ζ 1	Monod saturation	0.075	(0.05, 0.1)
ζ 2	Monod saturation	0.075	(0.05, 0.1)
c 1	Biofilm phage lysis	0.075	(0.05, 0.1)
c 2	Planktonic phage lysis	0.075	(0.05, 0.1)
p	Phage detachment rate	6.3	(2–8)
q	Phage reattachment rate	6.3	(2–8)
τ	Average latency time	6.3	(2–8)
b	Burst size	6.3	(2–8)
α	Infection decay rate	6.3	(2–8)

The range of values presented in this table is used for the parameter sensitivity analysis.

As illustrated in Figure 6, the PRCC indicates how variations in a specific parameter influence the output measure while accounting for the linear effects of other parameters. We calculated the PRCC values as Spearman (rank) partial correlations using the partialcorr function in MATLAB 2020, also determining their significance through uncorrelated p-values. The PRCC values range from − 1 to 1, with negative values suggesting an inverse relationship between the parameter and the outcome measure. Following the methodology of Marino et al. [28], we employed a z-test on the transformed PRCC values to rank the significant model parameters based on their relative sensitivity. Parameters with larger magnitudes were found to exert a stronger influence on the output measures.

Figure 6

PRCC. This is the partial rank correlation coefficient for the sensitivity analysis of the parameters. Source: This figure has been created by the authors.

Initially, we verify the monotonicity of the output measures as presented in Figures 7 and 8, which is confirmed for all parameters, with more emphasis on parameters that relate to bacteriophages and infected cells respectively. The PRCC analysis across these ranges yields consistent results. For the biofilm–phage model, we compute the PRCC for the following output measures: infected bacterial cells, bacteriophages in the planktonic phase, and 10% of the population of infected bacterial cells. Our findings show that β is not very sensitive to the population of infected cells because the infected class is more like a counter that captures how much of the bacterial cells are infected during the interaction, and this is not directly related to β but to ϕ 1 and ϕ 2 . Similarly, c 1 and c 2 are not sensitive to the bacteria in the planktonic phage, and this is the phage mortality rate in the biofilm. Since the bacteria and the phages are assumed to be region specific, these parameters do not affect the outputs of bacteria in the planktonic regions. Furthermore, ζ 1 and ζ 2 are not sensitive parameters to the infected population rather they are significantly sensitive to the bacteria population in the biofilm and planktonic regions respectively. For the parameters not showing monotonicity, we can attempt to solve the monotonicity by breaking down the graph into two monotonic regions, such that instead of the small range of outcome measures observed, and we would have considered truncating the range and looking at each truncated half separately. However, the current effect of these parameters on the overall outcome is not very significant.

Figure 7

Monotonicity check (infected population). This plot shows a monotonicity check of the parameters for the infected population. Source: This figure has been created by the authors.

Figure 8

Monotonicity check (bacteriophage population). This plot shows a monotonicity check of the parameters for the bacteriophages population in biofilm and planktonic regions. Source: This figure has been created by the authors.

5 Discussion

We have developed a bacteria–phage interaction model within a biofilm in a cystic fibrosis patient. The model considers bacteria in a biofilm and planktonic phase. The model assumes that the interactions are region specific, which means that the phages in the planktonic phase can only interact with the planktonic bacteria, while the phages in the biofilm can only interact with biofilm bacteria cells. The model speculates that the burst size could control the biofilm growth. In an effort to understand the state of the disease over time, we developed a stochastic model with which we could investigate the probabilities of reaching peak infection within a short time.

The model in this study can be easily adopted to investigate the effect of factors such as temperature and pH value on middle ear infection. Experimental study in previous studies [34,45] revealed that the pH of middle ear fluid collected from acute otitis media of children could affect biofilm formation, and biofilm formation is limited or completely absent under aerobic conditions, which is likely to happen; therefore, the current model in this study can be adopted with the inclusion of these specific factors to understand the interaction of phages and bacteria in middle ear infection.

Our assumptions and findings are consistent with the dynamics associated with biofilms. For instance, one of our main assumptions confirms that the phage interaction rate in the biofilm is different from the planktonic since the bacteria in the biofilm are dense and encased with extracellular polymeric substances, this is consistent with several in-vitro settings [2,4,5,12,18,23,44].

In connecting models to experiment, our model is not able to explain the biofilm occupancy, which will require the spatial components incorporated into the model, and the spatial structures can definitely be therapeutically relevant. For example, understanding the biofilm matrix and EPS will help to understand the actual interaction rates within the biofilm; the bacteria occupancy in the biofilm will help to determine if the phage interactions is at the biofilm surface, mimicking a lollipop-like degradation, or from within the biofilm, thus forming cavities. Another extension of this model are to (a) investigate a combination therapy that will involve antibiotics and immune response and (b) investigate the factors that influence biofilm formation and how they can be manipulated to prevent and eliminate biofilm-associated diseases in other areas of the body.

Acknowledgements

The authors appreciate the reviewers for their careful reading and valuable comments, which significantly enhanced the manuscript.

Funding information: The authors state no funding involved.
Author contributions: All authors contributed equally in analyzing the model, reviewing all the results, and preparing the manuscript. BOE conceptualized the study and designed the simulations. BOE and RNGRG performed the stability analysis.
Conflict of interest: The authors state no conflict of interest.
Ethical approval: The conducted research is not related to either human or animals use.
Data availability statement: This is a simulation paper. All model parameters needed to reproduce the study are given in the paper. No special software is required. The method used by us is described in detail in the literature as per references.

References

[1] Bañuelos, S., Gulbudak, H., Horn, M. A., Huang, Q., Nandi, A., Ryu, H., & Segal, R. (2020). Investigating the impact of combination phage and antibiotic therapy: A modeling study. Using Mathematics to Understand Biological Complexity. In: Association for Women in Mathematics Series, vol 22. Cham: Springer. DOI: https://doi.org/10.1007/978-3-030-57129-0_6. Search in Google Scholar

[2] Bardina, C., Spricigo, D. A., Cortés, P., & Llagostera, M. (2012). Significance of the bacteriophage treatment schedule in reducing salmonella colonization of poultry. Applied and Environmental Microbiology, 78(18), 6600–6607. DOI: https://doi.org/10.1128/aem.01257-12. Search in Google Scholar PubMed PubMed Central

[3] Beke, G., Stano, M., & Klucar, L. (2016). Modelling the interaction between bacteriophages and their bacterial hosts. Mathematical Biosciences, 279, 27–32. DOI: https://doi.org/10.1016/j.mbs.2016.06.009. Search in Google Scholar PubMed

[4] Bjarnsholt, T. (2013). The role of bacterial biofilms in chronic infections. APMIS, 121(s136), 1–58. DOI: http://dx.doi.org/10.1111/apm.12099. Search in Google Scholar PubMed

[5] Brüssow, H. (2013). Bacteriophage-host interaction: From splendid isolation into a messy reality. Current Opinion in Microbiology, 16(4), 500–506. DOI: https://doi.org/10.1016/j.mib.2013.04.007. Search in Google Scholar PubMed

[6] Carlton R. M. (1999). Phage therapy: past history and future prospects. Archivum immunologiae et therapiae experimentalis, 47(5), 267–274. Search in Google Scholar

[7] Cerca, N., Oliveira, R., & Azeredo, J. (2007). Susceptibility of staphylococcus epidermidis planktonic cells and biofilms to the lytic action of Staphylococcus bacteriophage K. Letters in Applied Microbiology, 45(3), 313–317. DOI: https://doi.org/10.1111/j.1472-765x.2007.02190.x. Search in Google Scholar

[8] D’Acunto, B., Frunzo, L., Klapper, I., Mattei, M. R., & Stoodley, P. (2019). Mathematical modeling of dispersal phenomenon in biofilms. Mathematical Biosciences, 307, 70–87. DOI: https://doi.org/10.1016/j.mbs.2018.07.009. Search in Google Scholar PubMed

[9] Delbrück, M. (1940). The growth of bacteriophage and lysis of the host. Journal of General Physiology, 23(5), 643–660. DOI: https://doi.org/10.1085/jgp.23.5.643. Search in Google Scholar PubMed PubMed Central

[10] Dieckmann, U. (2002). Adaptive dynamics of pathogen-host interactions. In Dieckmann U., Metz J. A. J., Sabelis M. W., & Sigmund K. editors, Adaptive Dynamics of Infectious Diseases: In Pursuit of Virulence Management, (pp. 39–59). chapter, Cambridge University Press, Cambridge. 10.1017/CBO9780511525728.006Search in Google Scholar

[11] De Paepe, M., & Taddei, F. (2006). Viruses’ life history: Towards a mechanistic basis of a trade-off between survival and reproduction among phages. PLoS Biology, 4(7), e193. DOI: https://doi.org/10.1371/journal.pbio.0040193. Search in Google Scholar PubMed PubMed Central

[12] Dkabrowska, K., & Abedon, S. T. (2019). Pharmacologically aware phage therapy: Pharmacodynamic and pharmacokinetic obstacles to phage antibacterial action in animal and human bodies. Microbiology and Molecular Biology Reviews, 83(4), 10–128. DOI: https://doi.org/10.1128/mmbr.00012-19. Search in Google Scholar

[13] Donlan, R. M. (2009). Preventing biofilms of clinically relevant organisms using bacteriophage. Trends in Microbiology, 17(2), 66–72. DOI: https://doi.org/10.1016/j.tim.2008.11.002. Search in Google Scholar PubMed

[14] Ellis, E. L., & Delbrück, M. (1939). The growth of bacteriophage. Journal of General Physiology, 22(3), 365–384. DOI: https://doi.org/10.1085/jgp.22.3.365. Search in Google Scholar PubMed PubMed Central

[15] Emerenini, B. O., & Eberl, H. J. (2022). Reactor scale modeling of quorum sensing induced biofilm dispersal. Applied Mathematics and Computation, 418, 126792. DOI: https://doi.org/10.1016/j.amc.2021.126792. Search in Google Scholar

[16] Emerenini, B. O., Hense, B. A., Kuttler, C., & Eberl, H. J. (2015). A mathematical model of quorum sensing induced Biofilm Detachment. PLoS One, 10(7), e0132385. DOI: https://doi.org/10.1371/journal.pone.0132385. Search in Google Scholar PubMed PubMed Central

[17] Emerenini, B. O., Ijioma, E. R., Wurscher, K., Reyes Grimaldo, R. N. G., & Williams, R. (2021). Mathematical modeling and analysis of influenza in-host infection dynamics. Letters in Biomathematics, 8(1), 229–253. DOI: https://doi.org/10.30707/lib8.1.1647878866.124006. Search in Google Scholar

[18] Eriksen, R. S., Mitarai, N., & Sneppen, K. (2020). Sustainability of spatially distributed bacteria–phage systems. Scientific Reports, 10(1), 3154. DOI: https://doi.org/10.1038/s41598-020-59635-7. Search in Google Scholar PubMed PubMed Central

[19] Ewald, J., Sieber, P., Garde, R., Lang, S. N., Schuster, S., & Ibrahim, B. (2019). Trends in mathematical modeling of host-Pathogen Interactions. Cellular and Molecular Life Sciences, 77(3), 467–480. DOI: https://doi.org/10.1007/s00018-019-03382-0. Search in Google Scholar PubMed PubMed Central

[20] Hadas, H., Einav, M., Fishov, I., & Zaritsky, A. (1997). Bacteriophage T4 development depends on the physiology of its host escherichia coli. Microbiology, 143(1), 179–185. DOI: https://doi.org/10.1099/00221287-143-1-179. Search in Google Scholar PubMed

[21] Hanlon, G. W., Denyer, S. P., Olliff, C. J., & Ibrahim, L. J. (2001). Reduction in exopolysaccharide viscosity as an aid to bacteriophage penetration through pseudomonas aeruginosa biofilms. Applied and Environmental Microbiology, 67(6), 2746–2753. DOI: https://doi.org/10.1128/aem.67.6.2746-2753.2001. Search in Google Scholar

[22] Hosseinidoust, Z., Tufenkji, N., & van de Ven, T. G. M. (2013). Formation of biofilms under phage predation: Considerations concerning a biofilm increase. Biofouling, 29(4), 457–468. DOI: https://doi.org/10.1080/08927014.2013.779370. Search in Google Scholar PubMed

[23] Hughes, K. A., Sutherland, I. W., & Jones, M. V. (1998). Biofilm susceptibility to bacteriophage attack: The role of phage-borne polysaccharide depolymerase. Microbiology, 144(11), 3039–3047. DOI: https://doi.org/10.1099/00221287-144-11-3039. Search in Google Scholar PubMed

[24] Jamal, M., Hussain, T., Das, C. R., & Andleeb, S. (2015). Characterization of siphoviridae phage Z and studying its efficacy against multidrug-resistant Klebsiella pneumoniae planktonic cells and biofilm. Journal of Medical Microbiology, 64(4), 454–462. DOI: https://doi.org/10.1099/jmm.0.000040. Search in Google Scholar PubMed

[25] Keeling, M. J., & Rohani, P. (2011). Modeling infectious diseases in humans and animals. Princeton University Press, Princeton. 10.2307/j.ctvcm4gk0Search in Google Scholar

[26] Kermack, W. O., & McKendrick, A. G. (1927). A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 115(772), 700–721. DOI: https://doi.org/10.1098/rspa.1927.0118. Search in Google Scholar

[27] Loc-Carrillo, C., & Abedon, S. T. (2011). Pros and cons of phage therapy. Bacteriophage, 1(2), 111–114. DOI: https://doi.org/10.4161/bact.1.2.14590. Search in Google Scholar PubMed PubMed Central

[28] Marino, S., Hogue, I. B., Ray, C. J., & Kirschner, D. E. (2008). A methodology for performing global uncertainty and sensitivity analysis in systems biology. Journal of Theoretical Biology, 254(1), 178–196. DOI: http://dx.doi.org/10.1016/j.jtbi.2008.04.011. Search in Google Scholar PubMed PubMed Central

[29] May, R. M., & Anderson, R. M. (1992). Infectious Diseases of Humans: Dynamics and Control. Oxford University Press, Oxford, UK. Search in Google Scholar

[30] McKay, M. D. (1992). Latin hypercube sampling as a tool in uncertainty analysis of computer models. Proceedings of the 24th Conference on Winter Simulation – WSC ’92 (pp. 557–564). DOI: https://doi.org/10.1145/167293.167637. Search in Google Scholar

[31] Mckay, M. D., Beckman, R. J., & Conover, W. J. (2000). A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 42(1), 55–61. DOI: http://dx.doi.org/10.1080/00401706.2000.10485979. Search in Google Scholar

[32] Musk Jr., D., & Hergenrother, P. (2006). Chemical countermeasures for the control of bacterial biofilms: Effective compounds and promising targets. Current Medicinal Chemistry, 13(18), 2163–2177. DOI: https://doi.org/10.2174/092986706777935212. Search in Google Scholar PubMed

[33] Ng, R. N., Tai, A. S., Chang, B. J., Stick, S. M., & Kicic, A. (2021). Overcoming challenges to make bacteriophage therapy standard clinical treatment practice for cystic fibrosis. Frontiers in Microbiology, 11, 593988. DOI: https://doi.org/10.3389/fmicb.2020.593988. Search in Google Scholar PubMed PubMed Central

[34] Osgood, R., Salamone, F., Diaz, A., Casey, J. R., Bajorski, P., & Pichichero, M. E. (2015). Effect of Ph and oxygen on biofilm formation in acute otitis media associated NTHi clinical isolates. The Laryngoscope, 125(9), 2204–2208. DOI: http://dx.doi.org/10.1002/lary.25102. Search in Google Scholar PubMed

[35] Parasion, S., Kwiatek, M., Gryko, R., Mizak, L., & Malm, A. (2014). Bacteriophages as an alternative strategy for fighting biofilm development. Polish Journal of Microbiology, 63(2), 137–145. 10.33073/pjm-2014-019Search in Google Scholar

[36] Potera, C. (2013). Phage renaissance: New hope against antibiotic resistance. Environmental Health Perspectives, 121(2), a48–a53. DOI: https://doi.org/10.1289/ehp.121-a48. Search in Google Scholar PubMed PubMed Central

[37] Ross, R. (1911). The Prevention of Malaria. John Murray, London. Search in Google Scholar

[38] Siettos, C. I., & Russo, L. (2013). Mathematical modeling of infectious disease dynamics. Virulence, 4(4), 295–306. DOI: https://doi.org/10.4161/viru.24041. Search in Google Scholar PubMed PubMed Central

[39] Sillankorva, S., Oliveira, R., Vieira, M. J., Sutherland, I., & Azeredo, J. (2004). Bacteriophage Φ S1 infection of pseudomonas fluorescens planktonic cells versus biofilms. Biofouling, 20(3), 133–138. DOI: https://doi.org/10.1080/08927010410001723834. Search in Google Scholar PubMed

[40] Sinha, S., Grewal, R. K., & Roy, S. (2020). Modeling phage–bacteria dynamics. Methods in Molecular Biology, 2131, 309–327. DOI: https://doi.org/10.1007/978-1-0716-0389-5_18. Search in Google Scholar PubMed

[41] van den Driessche P. (2017). Reproduction numbers of infectious disease models. Infectious Disease Modelling, 2(3), 288–303. DOI: https://doi.org/10.1016/j.idm.2017.06.002. Search in Google Scholar PubMed PubMed Central

[42] van den Driessche, P., & Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences, 180(1–2), 29–48. DOI: http://dx.doi.org/10.1016/s0025-5564(02)00108-6. Search in Google Scholar PubMed

[43] Wang, I.-N. (2006). Lysis timing and bacteriophage fitness. Genetics, 172(1), 17–26. DOI: https://doi.org/10.1534/genetics.105.045922. Search in Google Scholar PubMed PubMed Central

[44] Waters, E. M., Neill, D. R., Kaman, B., Sahota, J. S., Clokie, M. R., Winstanley, C., & Kadioglu, A. (2017). Phage therapy is highly effective against chronic lung infections with Pseudomonas aeruginosa. Thorax, 72(7), 666–667. DOI: http://dx.doi.org/10.1136/thoraxjnl-2016-209265. Search in Google Scholar PubMed PubMed Central

[45] Wright, A., Hawkins, C. H., Änggård, E. E., & Harper, D. R. (2009). A controlled clinical trial of a therapeutic bacteriophage preparation in chronic otitis due to antibiotic-resistant Pseudomonas aeruginosa; a preliminary report of efficacy. Clinical Otolaryngology, 34(4), 349–357. DOI: http://dx.doi.org/10.1111/j.1749-4486.2009.01973.x. Search in Google Scholar PubMed

Received: 2023-07-30

Revised: 2024-10-26

Accepted: 2024-12-04

Published Online: 2025-01-31

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/cmb-2024-0018

Keywords for this article

biofilm; phages; mathematical modeling; stochastic modeling

Creative Commons

BY 4.0