Mathematical modelling of CRISPR-Cas system effects on biofilm formation

ABSTRACT

Clustered regularly interspaced short palindromic repeats (CRISPR), linked with CRISPR associated (Cas) genes, can confer adaptive immunity to bacteria, against bacteriophage infections. Thus from a therapeutic standpoint, CRISPR immunity increases biofilm resistance to phage therapy. Recently, however, CRISPR-Cas genes have been implicated in reducing biofilm formation in lysogenized cells. Thus CRISPR immunity can have complex effects on phage–host–lysogen interactions, particularly in a biofilm. In this contribution, we develop and analyse a series of dynamical systems to elucidate and disentangle these interactions. Two competition models are used to study the effects of lysogens (first model) and CRISPR-immune bacteria (second model) in the biofilm. In the third model, the effect of delivering lysogens to a CRISPR-immune biofilm is investigated. Using standard analyses of equilibria, stability and bifurcations, our models predict that lysogens may be able to displace CRISPR-immune bacteria in a biofilm, and thus suggest strategies to eliminate phage-resistant biofilms.

Keywords:

• Bacteria
• bacteriophage
• biofilm
• CRISPR-Cas system
• lysogens

Table

Nomenclature
Parameters	Description	Dimensions
H	Density of wild-type (non-lysogenic) bacteria	cells cm^−2
$H_{\mathrm{L}}$	Density of lysogens	cells cm^−2
$C_{\mathrm{S}}$	Density of bacteria with CRISPR-immunity/spacer	cells cm^−2
$C_{\mathrm{L}}$	Density of lysogens with CRISPR system	cells cm^−2
V	Density of phage	phage cm^−2
t	Time	hr
K	Carrying capacity of bacteria in biofilm	cells cm^−2
α	Prophage induction	hr^−1
r	Bacterial growth rate	hr^−1
b	Burst size	phage cell^−1
β	Adsorption rate	cm^2 phage^−1 hr^−1
γ	Phage loss rate	hr^−1
θ	Adhesion rate per cell	hr^−1
η	Slouging off rate of non-lysogens	hr^−1
${\eta }_{\mathrm{L}}$	Sloughing off rate of lysogens	hr^−1
φ	Planktonic bacteria forming biofilm	cells cm^−2 hr^−1
$p_{\mathrm{L}}$	Probability of lysogenization	−
$p_{\mathrm{F}}$	Probability of CRISPR failure	−
$p_{\mathrm{D}}$	Probability of cell death	−

1. Introduction

The co-existence of bacteria and bacteriophage has been of prolonged interest to evolutionary biologists [1]. Temperate phages are of particular interest as these viruses can reproduce either through the lytic cycle, causing bacterial cell death, or the lysogenic cycle, allowing both phage and host to survive [2, 7, 18, 57]. An important factor in this co-existence is a mechanism known as clustered regularly interspaced short palindromic repeats (CRISPR) that not only provides adaptive immunity against viral infections but also helps to protect the bacterial genome from foreign mobile elements such as phage and plasmids [11, 51].

The building blocks of CRISPR systems were first identified 30 years ago, when inter-spaced DNA repeats were found in specific regions of the Escherichia coli genome [30]. Research efforts quantified the variation in these sequences, their lengths and positions, in bacterial and archaea genomes [5, 9, 39] before the name CRISPR was first applied [31]. CRISPR systems have now been identified in nearly 45% of bacterial strains and have been classified into three types that are further divided according to the phylogeny of the CRISPR sequences and their CRISPR-associated (Cas) genes [11]. Moreover, multiple CRISPR-Cas systems exist in several bacterial strains demonstrating that different types are compatible with each other within a single cell [16, 58].

The CRISPR adaptive defense mechanism proceeds in three main steps [11]. Firstly, the CRISPR-Cas system requires a specific sequence in phage DNA known as the protospacer that is acquired as a spacer into a designated position in the bacterial genome known as the CRISPR-locus, which lies next to the Cas genes [8, 22]. This acquisition is made possible by means of a short conserved sequence, present in the vicinity of protospacers, known as the protospacer adjacent motif (PAM) [17, 42]. CRISPR-loci have the capacity to store hundreds of spacers in the form of an array, while each spacer in the array is surrounded by short palindromic repeats known as CRISPR-repeats [34]. In the second step, the CRISPR-locus is transcribed and expressed as a single long mRNA that is cleaved into a single spacer and partial repeat known as CRISPR-RNA (crRNA) [24]. Depending upon the type of CRISPR system [11, 47], crRNA is associated with a Cas protein which provides a strong response against any RNA/DNA sequences matching the protospacer and begins to cleave these viral sequences with the help of Cas enzymes; this third step is known as interference. In the case of subsequent infections by the same type of phage, the crRNA associated with the Cas protein promptly responds and cleaves the phage genome. However, the presence of a limited number of crRNA/Cas complexes inside the cell can lead to CRISPR failure in the case of multiple simultaneous viral infections. Moreover, small variations in PAM can also lead to immune failure resulting in phage infection. In addition, rapidly mutating phage can escape this bacterial immunity by varying the protospacer or PAM sequences. A balance between phage diversification and CRISPR immunity can result in a stable co-existing community [4, 27, 54].

Biofilm formation is another protective mechanism used by bacterial colonies. For example, Pseudomonas aeruginosa preferentially exhibits an anaerobic biofilm mode of growth at human body temperature [55, 59]. In particular, when this pathogen colonizes wound infections (e.g. suppurative or purulent bacterial infections) or the lungs of cystic fibrosis patients [15, 49], it forms micro-colony structures which exhibit extremely high levels of antibiotic resistance [15]. Bacterial cells undergo profound phenotypic changes during the transition from free-floating planktonic to biofilm-associated cells [44].

New advancements in phage therapies are considered to be a promising way to eradicate some otherwise untreatable bacterial infections [41, 46]. Although the therapeutic effect of a virulent phage is typically more effective than a temperate phage [50], virulent phage therapy is not always possible and therefore the use of temperate phage is sometimes inevitable [12, 45, 48]. In particular, nearly 100% lysogenization can be achieved for P. aeruginosa when exposed to some viruses, [36], but these lysogens are normally immune to super-infection and plaque formation by the same phage [56].

Recent studies have revealed an effect of the type-I CRISPR-Cas system on the regulation of group behaviours in one strain of P. aeruginosa, PA14 [29, 60]. In particular, it has been demonstrated that PA14 is unable to form biofilm and loses swarming motility when CRISPR-Cas identifies a specific protospacer and PAM in a prophage sequence in the bacterial genome [29]. This result demonstrates that the CRISPR-Cas system, in addition to its role as an adaptive immune response, can regulate the genomic content of bacteria and help protect bacterial colonies from lysogenization; in fact, very few prophage are observed in P. aeruginosa genomes that carry CRISPR systems [19, 37]. The regulation of bacterial genomic content by means of the CRISPR-Cas system has been found in other bacterial strains as well, in which not only is lysogeny prevented, but also existing prophage is targeted, typically resulting in cell death due to genomic breakdown [19]. This self-destruction of lysogens may protect colonies by reducing the chance of prophage induction.

The effects of CRISPR-Cas systems on bacteria–phage interactions have been studied theoretically using a range of mathematical modelling approaches [14, 23, 26, 28, 35]. To date, these models have been developed for virulent phage assuming that CRISPR systems function as adaptive immune systems, without affecting other cellular processes such as biofilm formation [29, 60] or the cleavage of prophage [19]. In parallel with these efforts, a number of models have been developed to study the ecological effects of bacterial group behaviours, particularly biofilm formation [6, 20, 21, 32, 40]. Modelling efforts directed towards phage infection in biofilms are relatively rare [1, 52].

The goal of this work is to investigate the effects of a temperate phage infection in a bacterial biofilm. Our approach is to develop a series of models that help to isolate and disentangle the roles of bacterial hosts, either with or without a CRISPR system, lysogens, and phage. Not surprisingly, our models predict that CRISPR-immune bacteria in a biofilm, such as P. aeruginosa in cystic fibrosis, will be difficult to eradicate by phage therapy. However our results suggest that lysogens without a functioning CRISPR system may in some parameter regimes be able to invade and dominate the biofilm, offering a possible avenue for therapy.

2. Model

2.1. Model formulation

We develop a series of three models to investigate the effect of infection by temperate phage on a bacterial biofilm, for a host species either with or without a CRISPR-Cas system. Motivated by predominantly single species biofilms as seen in P. aeruginosa [29], all of the models consider cells of a single bacterial species within a biofilm. Each model has two bacterial populations and a population of bacteriophage as described in Figure 1. Although biofilms have complex and time-varying spatial structures [40, 43], we model each of these populations in aggregate, allowing us to focus on the ‘big picture’ questions of potential equilibria and stability. In the first model, the population dynamics of host cells, H, lysogenized host cells, $H_{\mathrm{L}}$ and phage, V, are studied in the absence of a CRISPR system. In the second model, a CRISPR system is introduced; we consider CRISPR-immune host cells, $C_{\mathrm{S}}$, and lysogens with a CRISPR system (but no phage-specific spacer) $C_{\mathrm{L}}$ (Model 2). In the third model, CRISPR-immune bacteria $C_{\mathrm{S}}$ are considered along with lysogens without a CRISPR system. The first case of Model 3 describes the biofilm in isolation whereas in the second case a constant population of planktonic (free) lysogens ${\hat{H}}_{\mathrm{L}}$ in the environment may join the biofilm. The common parameters of all three models are described as follows.

Figure 1. Bacterial populations considered in the models. H represents the population of wild-type bacteria, $H_{\mathrm{L}}$ represents the population of bacteria with prophage (lysogens), $C_{\mathrm{S}}$ is the population of bacteria with a CRISPR-Cas system, $C_{\mathrm{L}}$ are lysogens with a CRISPR-Cas system. The bacterial populations H and $H_{\mathrm{L}}$ are considered in Model 1, $C_S$ and $C_{\mathrm{L}}$ in Model 2 and $H_{\mathrm{L}}$ and $C_S$ in Model 3.

The bacterial populations are modelled as cell densities per unit area of biofilm substratum, cells/cm². These populations can increase logistically with a maximum growth rate r, but are limited by a fixed number of potential attachment sites, given by carrying capacity K cells/cm². This capacity corresponds to the total number of cells that could be supported in the biofilm per unit area of substratum. Although bacterial detachment from a biofilm is a complex spatial process, here we follow the Freter model [33] and assume that loss is simply linear in population size. This assumption might be appropriate, for example, when flow conditions are the dominant factor in detachment. Thus, bacteria detach from the biofilm with linear sloughing off rate η, which may differ for CRISPR and non-CRISPR bacteria as explained further in the sections to follow. Unlike η, because the bacterial populations are all the same species and differ only in the presence/absence of a CRISPR system or a single prophage sequence, r and K are common to all bacterial populations.

We assume that the virus is able to diffuse fairly freely through biofilm channels, yielding mass action attachment kinetics. Thus, $\beta HV$ gives the number of adsorption events per unit time, and infected bacterial cells typically undergo lysis, producing b daughter phage. The increase in the bacteriophage population is regulated by the rate of phage dissolution, γ; it is assumed that this phage loss includes phage particles diffusing from the biofilm. We also impose the standard assumption that the loss of phage via productive infection, $\beta HV$ is negligible compared to the overall clearance rate $\gamma V$. Using these parameters, three systems of ordinary differential equations (ODEs) are analysed below to predict which populations will persist at equilibrium for realistic parameter regimes.

2.2. Model 1

The first model is comparatively simple and based on the frequently studied phage–bacteria interaction when there is no CRISPR-Cas system present in the bacterial cells [3, 34, 36, 52]. Two bacterial populations are considered, the wild-type H and lysogens $H_{\mathrm{L}}$, along with a phage population V. The phage can attach to wild-type cells with adsorption rate β, affecting the host cell in one of two possible ways: (1) host cells H are lysogenized with probability $p_{\mathrm{L}}$, giving rise to lysogens $H_{\mathrm{L}}$ or (2) host cells are lysed with probability $(1-p_{\mathrm{L}})$ which results in the production of a burst of b phage.

In lysogens, prophage induction is possible at rate α which again results in cell lysis with burst size b. We assume lysogenization confers complete immunity to the virus, that is, $H_{\mathrm{L}}$ cannot be infected. The model is represented by the following set of ODEs: (1) $\begin{array}{rl}\frac{\mathrm{d}H\left(t\right)}{\mathrm{d}t}& =r\left(1-\frac{H+H_{\mathrm{L}}}{K}\right)H-\beta \mathrm{HV}-\eta H,\\ \frac{\mathrm{d}{H}_{\mathrm{L}}\left(t\right)}{\mathrm{d}t}& =r\left(1-\frac{H+H_{\mathrm{L}}}{K}\right)H_{\mathrm{L}}+p_{\mathrm{L}}\beta \mathrm{HV}-\alpha H_{\mathrm{L}}-\eta H_{\mathrm{L}},\\ \frac{\mathrm{d}V\left(t\right)}{\mathrm{d}t}& =b(1-p_{\mathrm{L}})\beta \mathrm{HV}-\gamma V+b\alpha H_{\mathrm{L}}.\end{array}$(1)

The purpose of this model is to study criteria for the stable existence of lysogens at equilibrium. This model clearly shows four steady states of which the first is the trivial equilibrium (TE${}_1$), in which all populations are zero. This is possible when the sloughing off rate η is sufficiently high to eradicate all bacteria from the biofilm. The second equilibrium is the phage free equilibrium (PFE${}_1$) in which only wild-type bacteria survive given by (2) ${\left[\begin{array}{c}H\\ H_{\mathrm{L}}\\ V\end{array}\right]}_{{\mathrm{PFE}}_1}=\left[\begin{array}{c}K(r-\eta )\\ 0\\ 0\end{array}\right].$(2)

The existence criterion for steady state PFE${}_1$ is $\eta <r$. The third equilibrium is the lysogenic equilibrium (LE${}_1$) in which lysogens survive but the wild-type bacteria die out. These lysogens produce prophage by induction which results in the survival of the phage population. However, these phage have no effect on the bacterial population since lysogens cannot be infected. The equilibrium LE${}_1$ is given by (3) ${\left[\begin{array}{c}H\\ H_{\mathrm{L}}\\ V\end{array}\right]}_{{\mathrm{LE}}_1}=\left[\begin{array}{c}0\\ K\frac{r-(\eta +\alpha )}{r}\\ \frac{R_{\alpha }}{r}(r-(\eta +\alpha \left)\right)\end{array}\right],$(3) where $R_{\alpha }=Kb\alpha /\gamma$. The existence condition for LE${}_1$ is $\eta <r-\alpha$ (note that PFE${}_1$ also exists when this condition holds). Finally, there is an all existing equilibrium (AEE${}_1$) in which all three populations are non-zero, given by, (4) ${\left[\begin{array}{c}H\\ H_{\mathrm{L}}\\ V\end{array}\right]}_{{\mathrm{AEE}}_1}=\left[\begin{array}{c}\frac{\overline{V}K{p}_{\mathrm{L}}}{R_{\beta }(\alpha -\overline{V})(1-p_{\mathrm{L}}))}\\ \frac{K}{R_{\beta }(1-p_{\mathrm{L}})}\\ \frac{K\overline{V}}{\beta }\end{array}\right],$(4) where $\overline{V}=(r-\eta )-r/R_{\beta }$ and $R_{\beta }=Kb\beta /\gamma$. The existence condition for AEE${}_1$ is $r\left(1-\frac{1}{R_{\beta }}\right)-\alpha <\eta <r\left(1-\frac{1}{R_{\beta }}\right).$

The eigenvalues of the Jacobian of each equilibrium provide stability conditions: these are provided in detail in Appendix S1. We find that TE${}_1$ is only stable when $r<\eta$, that is, if the sloughing off rate exceeds the maximum growth rate of the biofilm cells. For the local stability of PFE${}_1$, the existence condition $\eta <r$ ensures that the first eigenvalue is negative. Considering the second eigenvalue, we find the stability condition $\eta >r(1-1/R_{\beta })$ for PFE${}_1$, which holds when $R_{\beta }<1$. $R_{\beta }$ is analogous to a basic reproductive ratio for the phage, and this result implies that when $R_{\beta }<1$, the disease-free state PFE${}_1$ is stable.

The stability conditions for LE${}_1$ are closely related to the conditions for PFE${}_1$ with an additional dependence on the phage induction rate α; the condition for LE${}_1$ stability is $\eta <r(1-1/R_{\beta })-\alpha$. If we consider α to be sufficiently small, then stability switches between LE${}_1$ and PFE${}_1$ at $\eta =r(1-1/R_{\beta })$. Moreover, we see that $R_{\beta }>1$ must hold for the stability of LE${}_1$. The eigenvalues of the Jacobian evaluated at AEE${}_1$ are complicated analytical expressions. Therefore, the stability will be explored numerically in the next section. The existence and stability criteria for each steady state in Model 1 are summarized in Table 1.

Table 1. Existence and stability conditions for the steady states in Model 1.

Steady State	Existence Criteria	Stability Criteria
TE${}_1$	–	$\eta >r$
PFE${}_1$	$\eta <r$	$\eta >r\left(1-{\displaystyle \frac{1}{R_{\beta }}}\right)$
LE${}_1$	$\eta <r-\alpha$	$\eta <r\left(1-{\displaystyle \frac{1}{R_{\beta }}}\right)-\alpha$
AEE${}_1$	$r\left(1-{\displaystyle \frac{1}{R_{\beta }}}\right)-\alpha <\eta <r\left(1-{\displaystyle \frac{1}{R_{\beta }}}\right)$	Evaluated numerically

Note: The abbreviations used in this table are $R_{\beta }=Kb\beta /\gamma$, TE${}_1$: Trivial equilibrium, PFE${}_1$: Phage free equilibrium, LE${}_1$: Lysogenic equilibrium, AEE${}_1$: all existing equilibrium.

We note that, for biologically realistic parameter regimes, the phage induction rate α is likely to be very small relative to r and η. Thus, AEE${}_1$ rarely exists, the existence criteria for PFE${}_1$ and LE${}_1$ are effectively identical, and their stability criteria are complementary. The population can stably exist in three states: if $\eta >r$ the host cell population is not sustainable and the trivial equilibrium is stable. If $\eta <r$, the stability criterion for PFE${}_1$ determines whether the phage-free or lysogenic equilibrium is stable. In agreement with intuition, the PFE${}_1$ will be stable for $R_{\beta }<1$.

2.3. Model 2

In this model, we examine how the situation above might differ if the bacterial strains involved have CRISPR immunity. We thus assume that both the wild-type and lysogenized cells have a CRISPR-Cas system, and denote these populations $C_{\mathrm{S}}$ and $C_{\mathrm{L}}$, respectively. The first population is assumed to have previously acquired the protospacer of phage DNA; these $C_{\mathrm{S}}$ cells are therefore CRISPR-immune. CRISPR-immune bacteria have the ability to resist future infections by the same phage, however, there is a small possibility of failure. Therefore, in the case of phage adsorption by $C_{\mathrm{S}}$, the probability of CRISPR failure $p_{\mathrm{F}}$ leads the cell to lysis.

For lysogens with a CRISPR system, we note results from a recent experimental study [19] which demonstrates that bacteriophage infection of previously lysogenized bacteria may lead to cell death, if the CRISPR-Cas proteins acquire the protospacer of the infecting phage. In this case, the same adsorption rate (β) is considered for the interaction between the bacteriophage V and lysogens $C_{\mathrm{L}}$, however, this interaction is lethal to the cell with probability $p_{\mathrm{D}}$. Finally, lysogens have a CRISPR-Cas system which regulates their group behaviour [29] so that they lose the expression of polysaccharides essential for biofilm formation. Therefore, the sloughing off rate in $C_{\mathrm{L}}$ increases (${\eta }_{\mathrm{L}}>\eta$ ). This dynamical system can be written as (5) $\begin{array}{rl}\frac{\mathrm{d}{C}_{\mathrm{S}}\left(t\right)}{\mathrm{d}t}& =r\left(1-\frac{C_{\mathrm{S}}+C_{\mathrm{L}}}{K}\right)C_{\mathrm{S}}-p_{\mathrm{F}}\beta C_{\mathrm{S}}V-\eta C_{\mathrm{S}},\\ \frac{\mathrm{d}{C}_{\mathrm{L}}\left(t\right)}{\mathrm{d}t}& =r\left(1-\frac{C_{\mathrm{S}}+C_{\mathrm{L}}}{K}\right)C_{\mathrm{L}}-p_{\mathrm{D}}\beta C_{\mathrm{L}}V-\alpha C_{\mathrm{L}}-{\eta }_{\mathrm{L}}{C}_{\mathrm{L}},\\ \frac{\mathrm{d}V\left(t\right)}{\mathrm{d}t}& =b{p}_{\mathrm{F}}\beta C_{\mathrm{S}}V-\gamma V+b\alpha C_{\mathrm{L}}.\end{array}$(5)

We begin by examining the existence conditions for equilibria of this model. Apart from the trivial equilibrium (TE${}_2=(0,0,0)$), there are four more equilibrium states: the phage free equilibrium (PFE${}_2$) in which only $C_{\mathrm{S}}$ survives, the lysogenic equilibrium (LE${}_2$) in which $C_{\mathrm{L}}$ and V survive, the CRISPR equilibrium (CE${}_2$) in which $C_{\mathrm{S}}$ and V survive and the all existing equilibrium (AEE${}_2$) in which all three populations exist. The first three non-trivial equilibria can be written as (6) $\left[\begin{array}{c}{C}_{\mathrm{S}}\\ C_{\mathrm{L}}\\ V\end{array}\right]={\left[\begin{array}{c}\frac{K(r-\eta )}{r}\\ 0\\ 0\end{array}\right]}_{PF{E}_2},{\left[\begin{array}{c}0\\ K\frac{r-({\eta }_{\mathrm{L}}+\alpha )}{R_{\beta }\alpha p_{\mathrm{D}}+r}\\ R_{\alpha }\frac{r-({\eta }_{\mathrm{L}}+\alpha )}{R_{\beta }\alpha p_{\mathrm{D}}+r}\end{array}\right]}_{L{E}_2},{\left[\begin{array}{c}\frac{K}{p_{\mathrm{F}}{R}_{\beta }}\\ 0\\ \frac{(r-\eta )-\frac{r}{p_{\mathrm{F}}{R}_{\beta }}}{\beta p_{\mathrm{F}}}\end{array}\right]}_{C{E}_2}.$(6)

The existence conditions for the phage free and lysogenic equilibria are similar to those obtained from the previous model and are given in Table 2. The equilibrium state for AEE${}_2$ is given by (7) $\left[\begin{array}{c}{C}_{\mathrm{S}}\\ C_{\mathrm{L}}\\ V\end{array}\right]={\left[\begin{array}{c}K\frac{\alpha \left(p_{\mathrm{D}}\right(r-\eta )-p_{\mathrm{F}}(r-(\alpha +{\eta }_{\mathrm{L}})\left)\right)+{\displaystyle \frac{r}{R_{\beta }}}(\alpha +{\eta }_{\mathrm{L}}-\eta )}{r\left(p_{\mathrm{F}}\right({\eta }_{\mathrm{L}}-\eta )+\alpha p_{\mathrm{D}})}\\ K{p}_{\mathrm{F}}(\alpha +{\eta }_{\mathrm{L}}-\eta )\frac{p_{\mathrm{F}}\left(r\left(1-{\displaystyle \frac{1}{R_{\beta }{p}_{\mathrm{F}}}}\right)-(\alpha +{\eta }_{\mathrm{L}})\right)-p_{\mathrm{D}}\left(r\left(1-{\displaystyle \frac{1}{R_{\beta }{p}_{\mathrm{F}}}}\right)-\eta \right)}{r(p_{\mathrm{F}}-p_{\mathrm{D}})\left(p_{\mathrm{F}}\right({\eta }_{\mathrm{L}}-\eta )+\alpha p_{\mathrm{D}})}\\ \frac{\alpha +{\eta }_{\mathrm{L}}-\eta }{\beta (p_{\mathrm{F}}-p_{\mathrm{D}})}\end{array}\right]}_{{\mathrm{AEE}}_2}.$(7)

Table 2. Existence and stability conditions of the steady states in Model 2.

Steady state	Existence criteria	Stability criteria
TE${}_2$	–	$\eta >r$
PFE${}_2$	$\eta <r$	$\eta >r\left(1-{\displaystyle \frac{1}{R_{\beta }{p}_{\mathrm{F}}}}\right)$
LE${}_2$	$\eta <{\eta }_{\mathrm{L}}<r-\alpha$	$\eta >r-(r-(\alpha +{\eta }_{\mathrm{L}}\left)\right){\displaystyle \frac{R_{\beta }\alpha p_{\mathrm{F}}+r}{R_{\beta }\alpha p_{\mathrm{D}}+r}}$
CE${}_2$	$\eta <r\left(1-{\displaystyle \frac{1}{R_{\beta }{p}_{\mathrm{F}}}}\right)$	$\eta <r\left(1-{\displaystyle \frac{1}{R_{\beta }{p}_{\mathrm{F}}}}\right)\left(1-{\displaystyle \frac{p_{\mathrm{F}}}{p_{\mathrm{D}}}}\right)+{\displaystyle \frac{p_{\mathrm{F}}}{p_{\mathrm{D}}}}(\alpha +{\eta }_{\mathrm{L}})$
AEE${}_2$	$\eta <r-(r-(\alpha +{\eta }_{\mathrm{L}}\left)\right){\displaystyle \frac{R_{\beta }\alpha p_{\mathrm{F}}+r}{R_{\beta }\alpha p_{\mathrm{D}}+r}}$
	$\eta >r\left(1-{\displaystyle \frac{1}{R_{\beta }{p}_{\mathrm{F}}}}\right)\left(1-{\displaystyle \frac{p_{\mathrm{F}}}{p_{\mathrm{D}}}}\right)+{\displaystyle \frac{p_{\mathrm{F}}}{p_{\mathrm{D}}}}(\alpha +{\eta }_{\mathrm{L}})$, $p_{\mathrm{F}}>p_{\mathrm{D}}$	Evaluated numerically

Note: The abbreviations used in this table are $R_{\beta }=Kb\beta /\gamma$, TE${}_2$: Trivial equilibrium, PFE${}_2$: Phage free equilibrium, LE${}_2$: Lysogenic equilibrium, CE${}_2$: CRISPR equilibrium and AEE${}_2$: All existing equilibrium.

The existence conditions for the above equilibrium are quite complex as shown in Table 2, however, one necessary condition is $p_{\mathrm{F}}>p_{\mathrm{D}}$. Note that this is the only equilibrium whose existence depends on the relation between the probabilities of CRISPR failure and cell death.

All these conditions admit the possibility of co-existence of equilibria and we thus consider stability conditions from the eigenvalues of the Jacobian, as provided in Appendix S2. In brief, for PFE${}_2$ we find an analogous condition to that found in previous model for PFE${}_1$; in particular, PFE${}_2$ is locally stable under the condition $\eta >r(1-1/R_{\beta }{p}_{\mathrm{F}})$. For CE${}_2$, the second and third eigenvalues are complex conjugates with negative real parts whenever CE${}_2$ exists, and we find a necessary condition for the negativity of the first eigenvalue. Similarly, for LE${}_2$, the second and third eigenvalues are negative whenever the equilibrium exists but we derive a condition for the negativity of the first eigenvalue. Finally, the eigenvalues for the Jacobian matrix of AEE${}_2$ are sufficiently complicated expressions that we will use numerical methods to assess stability. Table 2 shows a summary of existence and stability conditions for Model 2.

We find that the conditions for the stable existence of the trivial and phage-free equilibria are identical to those described for Model 1, with $R_{\beta }$ replaced by $R_{\beta }{p}_{\mathrm{F}}$. Thus again if the sloughing off rate exceeds the growth rate, the biofilm cannot sustain itself. If the biofilm can sustain itself, we find that the phage population cannot be sustained if $R_{\beta }{p}_{\mathrm{F}}<1$. If both the bacterial cells and phage are sustainable, however, three equilibrium states are possible. Once again, though, we note that α is very small, and thus the lysogenic equilibrium is rarely stable, since $\eta <{\eta }_{\mathrm{L}}$ by definition. The all existing equilibrium again has a possibly narrow parameter range for existence, and thus the model predicts that for biologically relevant parameter values, the CRISPR equilibrium is most likely to be observed. This is illustrated further in the numerical work to follow.

2.4. Model 3

This model consists of a population of CRISPR-immune bacteria, $C_{\mathrm{S}}$, along with non-CRISPR lysogens, $H_{\mathrm{L}}$, which continually contribute to the phage population V in the biofilm via induction. The idea here is to study the impact of non-CRISPR lysogens on the CRISPR-immune bacteria. The populations and parameters are the same as described for the previous models, with the exception of the biofilm formation rate φ. This parameter represents an external source of planktonic lysogens in the environment, that may join the biofilm. Including φ allows us to investigate the possibility of using lysogens to deliver phage therapy, with the long-term goal of eradicating the CRISPR-immune bacteria.

The model assumes that planktonic bacteria are present in the environment around the biofilm, allowing planktonic bacteria $\left(H_{\mathrm{L}}\right)$ to form biofilm at a specific rate constant of adhesion while the medium carrying planktonic bacteria flows into and out of the system, including the sloughed off bacteria from the biofilm. We use φ to denote the maximum rate of biofilm formation by planktonic bacteria; this rate is reduced by the carrying capacity such that the net biofilm formation rate depends on the attachment sites available in the biofilm. We assume that bacteria are lost from the biofilm independent of the density of planktonic bacteria; this sloughing off rate is assumed to be the same for both $C_{\mathrm{S}}$ and $H_{\mathrm{L}}$. The model is represented by the following equations: (8) $\begin{array}{rl}\frac{\mathrm{d}{C}_{\mathrm{S}}\left(t\right)}{\mathrm{d}t}& =r\left(1-\frac{C_{\mathrm{S}}+H_{\mathrm{L}}}{K}\right)C_{\mathrm{S}}-p_{\mathrm{F}}\beta C_{\mathrm{S}}V-\eta C_{\mathrm{S},}\\ \frac{\mathrm{d}{H}_{\mathrm{L}}\left(t\right)}{\mathrm{d}t}& =\left(1-\frac{C_{\mathrm{S}}+H_{\mathrm{L}}}{K}\right)(r{H}_{\mathrm{L}}+\phi )-\alpha H_{\mathrm{L}}-\eta H_{\mathrm{L}},\\ \frac{\mathrm{d}V\left(t\right)}{\mathrm{d}t}& =b{p}_{\mathrm{F}}\beta C_{\mathrm{S}}V-\gamma V+b\alpha H_{\mathrm{L}}.\end{array}$(8)

The above system is studied in two cases: (1) when there are no planktonic lysogens in the environment ($\phi =0$) and (2) when there are lysogens in the environment that may join the biofilm ($\phi >0$).

Case (a): $\mathit{\phi }\mathbf{=}\mathbf{0}$

When the biofilm formation term $\phi =0$, this model becomes relatively simple, yielding four steady states, three of which are defined in the previous model, i.e. the phage-free equilibrium ($PF{E}_3$), trivial equilibrium ($T{E}_3$) and CRISPR equilibrium ($C{E}_3$) with the population of lysogens $H_{\mathrm{L}}$ in place of $C_{\mathrm{L}}$, while the fourth equilibrium, i.e. the lysogenic $L{E}_3$, is as defined for $L{E}_1$ with $C_{\mathrm{S}}=0$ instead of H=0.

The eigenvalues of the Jacobian evaluated at these equilibria are given in Appendix S3, and the existence and stability conditions are summarized in Table 3. There are similarities in the stability conditions of each equilibrium of this model with those described for the previous models. In particular the stability conditions for the trivial and phage-free equilibria are the same as previously described for Model 2. In this case since α is small, the CRISPR equilibrium has a narrow range of stability, and biofilms are thus most likely to exist either at the phage-free or lysogenic equilibria.

Table 3. Existence and stability conditions of the steady states in Model 3.

Steady state	Existence criteria	Stability criteria
TE${}_3$	–	$\eta >r$
PFE${}_3$	$\eta <r$	$\eta >r\left(1-{\displaystyle \frac{1}{R_{\beta }{p}_{\mathrm{F}}}}\right)$
LE${}_3$	$\eta <r-\alpha$	$\eta <r\left(1-{\displaystyle \frac{1}{R_{\beta }{p}_{\mathrm{F}}}}\right)-\alpha$
CE${}_3$	$\eta <r\left(1-{\displaystyle \frac{1}{R_{\beta }{p}_{\mathrm{F}}}}\right)$	$r\left(1-{\displaystyle \frac{1}{p_{\mathrm{F}}{R}_{\beta }}}\right)-\alpha <\eta <r\left(1-{\displaystyle \frac{1}{p_{\mathrm{F}}{R}_{\beta }}}\right)$
LE${}_{\phi }$	−	$\phi >{\displaystyle \frac{K\alpha (r-\eta )}{\eta }}$ (sufficient condition)
AEE${}_{\phi }$	$\phi >{\displaystyle \frac{\alpha \gamma K(r-\eta )(\gamma r-Kb\beta p_{\mathrm{F}}(r-\eta -\alpha \left)\right)}{(Kb\alpha \beta p_{\mathrm{F}}+\gamma r)(Kb\alpha \beta p_{\mathrm{F}}+\eta \gamma )}}$
	$\phi <K\alpha \left(1-{\displaystyle \frac{r}{\alpha +\eta }}\left(1-{\displaystyle \frac{\gamma }{Kb\beta p_{\mathrm{F}}}}\right)\right)$	Evaluated numerically

Notes: The abbreviations used in this table are $R_{\beta }=Kb\beta /\gamma$, TE${}_3$: Trivial equilibrium, PFE${}_3$: Phage free equilibrium, LE${}_3$ and LE${}_{\phi }$: Lysogenic equilibrium, CE${}_3$: CRISPR Equilibrium and AEE${}_3$ and AEE${}_{\phi }$: All existing equilibrium.

2.4.1. Case (b): $\mathit{\phi }\mathbf{>}\mathbf{0}$

This situation changes when planktonic lysogens in the environment are able to join the biofilm. When $\phi >0$, only two equilibria remain. One is the lysogenic equilibrium (LE${}_{\phi }$) and the other is the all existing equilibrium (AEE${}_{\phi }$). Equilibrium LE${}_{\phi }$ is given by (9) $\left[\begin{array}{c}{C}_{\mathrm{S}}\\ H_{\mathrm{L}}\\ V\end{array}\right]={\left[\begin{array}{c}0\\ \frac{1}{2r}\left(K\right(r-\eta -\alpha )-\phi +\sqrt{\left(K\right(r-\eta -\alpha )-\phi {)}^2+4Kr\phi ))}\\ \frac{R_{\beta }}{2Kr}\left(K\right(r-\eta -\alpha )-\phi +\sqrt{\left(K\right(r-\eta -\alpha )-\phi {)}^2+4Kr\phi ))}\end{array}\right]}_{L{E}_{\phi }}.$(9)

It can be observed that the lysogens and phage population are always positive because the discriminant is always positive and greater than the expression outside the square-root.

For the stability of LE${}_{\phi }$, a complex condition arises from the first eigenvalue, as given in Appendix S3. However, the following condition is sufficient to ensure the stability of LE${}_{\phi }$: $K\alpha (r-\eta )-\eta \phi <0$ which implies that $\phi >K\alpha (r-\eta )/\eta$. This condition is further analysed in the numerical section to follow, in which the parameter space is explored to delineate the region of stability.

The second equilibrium state AEE${}_{\phi }$ is given by (10) $\left[\begin{array}{c}{C}_{\mathrm{S}}\\ H_{\mathrm{L}}\\ V\end{array}\right]={\left[\begin{array}{c}\frac{Kb\alpha (r-\eta -D_4\beta p_{\mathrm{F}})-D_4\gamma r}{br(\alpha -D_4\beta p_{\mathrm{F}})}\\ \frac{\phi (\eta +\beta p_{\mathrm{F}}{D}_4)}{r(\alpha -D_4\beta p_{\mathrm{F}})}\\ D_4\end{array}\right]}_{{\mathrm{AEE}}_{\phi }},$(10) where $C_4=\gamma r-b\beta p_{\mathrm{F}}(K(r-\eta )+\phi )$ and $D_4=(-C_4+\sqrt{C_4^2+4K{b}^2{\beta }^2\eta p_{\mathrm{F}}^2\phi })/(2Kb{\beta }^2{p}_{\mathrm{F}}^2)$. Since $D_4$ is always positive, we have two possible existence conditions: $\begin{array}{rl}& \left(1\right)\alpha -D_4\beta p_{\mathrm{F}}>0⟹\phi <K\alpha \left(1-\frac{r}{\alpha +\eta }\left(1-\frac{R_{\beta }}{p_{\mathrm{F}}}\right)\right)\mathrm{and}\\ & \left(2\right)Kb\alpha (r-\eta -D_4\beta p_{\mathrm{F}})-D_4\gamma r>0⟹\phi >\frac{K\alpha (r-\eta )(r-R_{\beta }{p}_{\mathrm{F}}(r-\eta -\alpha \left)\right)}{(R_{\beta }\alpha p_{\mathrm{F}}+r)(R_{\beta }\alpha p_{\mathrm{F}}+\eta )}.\end{array}$

Once again the eigenvalues of the Jacobian evaluated at AEE${}_{\phi }$ are not compact expressions, and stability will be explored numerically in the next section.

Model 3 has been developed to investigate the possibility of using lysogenized bacteria as a delivery mechanism in phage therapy. Without the addition of lysogenized bacteria ($\phi =0$), the most likely stable equilibrium for the realistic parameter values is PFE${}_3$, particularly when $R_{\beta }<1/p_{\mathrm{F}}$; this corresponds to a biofilm composed of bacteria that are immune to the phage. However, the stability conditions for the two equilibria when $\phi >0$ imply that when the rate of biofilm formation by lysogenized bacteria is sufficiently high, the population of CRISPR-immune bacteria may be eliminated from the biofilm.

3. Numerical simulations

Numerical simulations for the above models have been performed to illustrate the impact of parameter values on the existence and stability of the steady states. In addition, the eigenvalues for the Jacobian at two all existing equilibria (AEE${}_1$ and AEE${}_{\phi }$) which were not provided analytically are illustrated numerically.

3.1. Parameter values

Baseline parameters were obtained through a review of experimental results in the literature as well as parameter values used in previous mathematical modelling studies (Table 1) [3, 14, 21, 33, 36, 38]. The rate of bacterial growth r=1 hr⁻¹ is the maximum growth rate [1, 14] in the logistic growth term; growth slows at higher densities. The carrying capacity of the bacterial population inside the biofilm is the total attachment sites available per unit area. We take $K=5\times {10}^6$ cells cm⁻² as chosen by Freter [21] and followed by Ballyk et al. [6] and Mašić and Eberl [40].

We assume lysogens leave the cell population by means of prophage induction at the rate $\alpha ={10}^{-5}$ hr⁻¹, producing b=200 phage copies per cell through lysis. Phage can infect bacteria at the rate $\beta ={10}^{-7}$ cm² phage⁻¹ hr⁻¹ and produce b=200 phage copies per cell with probability $1-p_{\mathrm{L}}$, where $p_{\mathrm{L}}=0.05$ is the probability of bacteria gaining prophage. The phage loss rate is assumed to be $\gamma =0.05$ hr⁻¹ which includes phage dissolution and sloughing off from the biofilm.

Biofilm formation and sloughing off are two important phenomena in the study of bacterial group behaviour and these are frequently discussed in the literature [6, 10, 21, 40]. In this article, the sloughing off rate of bacteria is considered to be the same in all models, except that lysogens in Model 2, $C_{\mathrm{L}}$, leave the biofilm more rapidly due to the CRISPR response. The baseline parameter value for the sloughing off rate for wild-type bacteria is taken to be $\eta =0.1$ hr⁻¹ and for lysogens with a CRISPR system is given by ${\eta }_{\mathrm{L}}=\sqrt{\eta }$ so that ${\eta }_{\mathrm{L}}>\eta$ when $0<\eta <1$. The parameter value for the biofilm formation rate φ is found from Freter's model by assuming that the density of bacteria present near the biofilm is constant as $K_2={10}^8$ cells cm⁻³. The adhesion rate given in Freter's model, i.e. $\tilde{\theta }={10}^{-5}\mathrm{l}{\mathrm{hr}}^{-1}{\mathrm{g}}^{-1}$, is converted to $\theta =\tilde{\theta }{K}_2$ hr⁻¹ $\approx 5.56\times {10}^{-7}{\mathrm{hr}}^{-1}$ by assuming that 1 gram of bacterial mass can contain $N=1.8\times {10}^{12}$ cells at maximum. Therefore, the population forming biofilm per hour is given by $\phi =\theta K\mathrm{cells}{\mathrm{cm}}^{-2}{\mathrm{hr}}^{-1}$, where K is the carrying capacity of the biofilm.

The CRISPR-Cas effects on the wild-type bacteria and lysogens are parametrized so that they reflect biologically plausible dynamics. Since CRISPR immunity is considered highly efficient with a small probability of failure, we take $p_{\mathrm{F}}={10}^{-4}$ [14]. The probability of failure is a key measure in the models, allowing us to study phage infection in CRISPR-Cas bacteria. Failure is possible since the diversity in CRISPR arrays alters the effectiveness of the CRISPR-Cas system. In particular, CRISPR spacers may become less effective over time and therefore older spacers are preferentially replaced [25, 53]. The probability of death $p_{\mathrm{D}}$ depends on the adsorption rate of lysogens in the presence of CRISPR-Cas, the presence of repressors for the phage and the diversity of Cas proteins inside CRISPR-immune bacteria. In this study, the same adsorption rate is assumed for lysogens and non-lysogens, while $p_{\mathrm{D}}=0.1$ is used to model CRISPR-Cas initiated cell death (Table 4).

Table 4. Baseline values for the parameters used in this study.

Parameters	values	Source
K	$5\times {10}^6$ cells cm^−2	[6, 21, 40]
α	10^−5 hr^−1	[13, 36]
r	1 hr^−1	[1]
b	200 phage cell^−1	[1]
β	10^−7 cm^2 phage^−1 hr^−1	[1]
γ	0.05 hr^−1	[1]
η	0.1 hr^−1	[21]
${\eta }_{\mathrm{L}}$	$\sqrt{\eta }$ hr^−1	–
φ	2.78 cells cm^−2 hr^−1	see text
$p_{\mathrm{F}}$	10^−4	[14]
$p_{\mathrm{L}}$	0.05	–
$p_{\mathrm{D}}$	0.1	–

3.2. Baseline analysis

In Figure 2, the population dynamics of each of these models are illustrated for the same baseline parameter values and initial conditions. The simulations are run over long times to illustrate stable equilibrium states, and in each case the numerical results were validated by comparison with the analytical expressions. We plot both population densities and time on a log axis, such that initial transients are visible.

Figure 2. Dynamics of phage and bacterial populations inside the biofilm at baseline parametric values. H and $C_{\mathrm{S}}$ represent populations of wild-type bacteria and CRISPR-immune bacteria (dashed lines), $C_{\mathrm{L}}$ and $H_{\mathrm{L}}$ represent lysogens with and without CRISPR systems (dots) and V is the population of bacteriophage (solid lines). Model 1 approaches the lysogenic equilibrium LE${}_1$, Model 2 and Model 3 (case (a)) approach the phage-free equilibria (PFE${}_2$ and PFE${}_3$) whereas Model 3 (case (b)) approaches the all existing equilibrium (AEE${}_{\phi }$).

In the first model, the baseline scenario produces a stable lysogenic equilibrium (without CRISPR), in which the phage population is maintained in the biofilm through induction. This is a classical model which supports the idea of co-existence of temperate phage and their bacterial hosts. This result should generalize to the case of many phage types, since bacteria are capable of carrying several prophage as part of their genome.

In the second model, the CRISPR bacteria $C_{\mathrm{S}}$ clearly dominate, since CRISPR immunity makes these cells resistant to the phage. The CRISPR effects on the lysogens $C_{\mathrm{L}}$, i.e. the death of lysogens or increase in the sloughing-off rate, make it possible to eliminate these infected bacteria from the biofilm; their loss is followed by the loss of the phage population. The second model demonstrates that CRISPR systems will disrupt the coexistence observed in Model 1, and predicts that for bacteria with CRISPR systems, stable biofilms are likely to be composed of CRISPR-immune bacteria only.

The populations in the first two models reach stable equilibria fairly quickly, whereas simulation results for the third model show that the populations take longer to reach their stable states. In the absence of an external source of lysogens (Model 3a), CRISPR-immune bacteria dominate, since the CRISPR system helps the bacteria to maintain the phage-free equilibrium after eliminating non-CRISPR lysogens from the biofilm. When the external source of lysogens is present (Model 3b), we observe stability of the all existing equilibrium for these baseline parameter values. This is because a high biofilm formation rate is required to meet the sufficient condition for the stability of LE${}_{\phi }$, thus eliminating the CRISPR cells from the biofilm. To investigate this further, we turn to parametric analysis to explore scenarios in which lysogens are able to replace CRISPR-immune bacteria in the biofilm. The second goal of the parametric analysis is to find a suitable conditions at which, after the eradication of CRISPR bacteria from the biofilm, the lysogens themselves could be eradicated.

3.3. Parametric analysis

In this section we illustrate several bifurcations that may occur in biologically realistic parameter regimes. We first illustrate the changes in equilibria and stability that result from variation in the adsorption rate, β. Variation in the prophage induction rate, α, phage loss rate, γ, and sloughing off rate η have similar qualitative effects and are shown in Appendix S4. Finally, we examine the interesting case of variation in the biofilm formation rate, φ.

3.3.1. Adsorption rate β

Variation in the adsorption rate constant affects the stability of equilibrium states by varying $R_{\beta }$. Figure 3 shows that LE${}_1$ remains stable for a wide range of β values in Model 1, while the other models each switch stability between two equilibrium states at high rates of adsorption. In Model 2, PFE${}_2$ is replaced by CE${}_2$ at high adsorption rates. This is due to the strong immune response of CRISPR bacteria, such that even for high adsorption rates CRISPR-immune bacteria survive, co-existing with the phage.

Figure 3. Population densities of bacteria and bacteriophage at stable (solid lines) and unstable (dashed lines) equilibrium states against adsorption rate constant β in the above models. Each colour represents a unique population while the line width is increased to visualize overlapping populations. Bacterial populations are H (green) and $H_{\mathrm{L}}$ (blue) in Model 1, $C_{\mathrm{S}}$ (black) and $C_{\mathrm{L}}$ (purple) in Model 2 and $C_{\mathrm{S}}$ and $H_{\mathrm{L}}$ in Model 3 whereas the phage population is represented by V (red) in all models.

In Model 3, LE${}_2$ is stable at high β, replacing PFE${}_3$ in case (a) and AEE${}_{\phi }$ in case (b). We note that the biofilm formation rate φ does not have a strong influence here, since almost same adsorption rate β is required in both cases to eradicate CRISPR-immune bacteria.

Qualitatively, these results demonstrate that for biofilms without CRISPR systems (Model 1), lysogens and phage are predicted to coexist over a wide parameter range. In contrast, for biofilms with CRISPR (Model 2), the stable equilibrium state is typically dominated by CRISPR-immune bacteria. If CRISPR-capable lysogens are replaced by lysogens without a functioning CRISPR system (Model 3, a and b), again CRISPR-immune bacteria typically dominate, although for parameter regimes that are extremely favourable to phage reproduction, the lysogenic equilibrium may be stable. As shown in Appendix S4, these qualitative conclusions also hold for variations in parameters α, η and γ.

3.3.2. Biofilm formation rate φ

In Model 3, the rate of biofilm formation by planktonic lysogens is an important parameter that can be regulated either by increasing the density of lysogens near the biofilm or by increasing the adhesion rate. We illustrate the effect of varying this parameter in panel (a) of Figure 4. As the value φ increases, the lysogen population increases along with the phage, while the population of CRISPR-immune bacteria $C_{\mathrm{S}}$ decreases in the biofilm, and is eventually eliminated at high values of φ. Though a relatively large formation rate is required to eliminate $C_{\mathrm{S}}$, this threshold value depends on other model parameters. In particular, the critical φ value to ensure the stability of the lysogenic equilibrium in Model 3 (case (b)) is give by $\phi >K\alpha \left((r-\eta )/\eta \right)$. In the lower three panels of Figure 4, we plot this threshold against α, η and the carrying capacity K. Clearly, an increase in the value of α and K require an increase in the value of φ to make the lysogenic equilibrium stable, whereas an increase in the sloughing off rate η reduces the threshold value of φ.

Figure 4. (a) Population densities of bacteria and bacteriophage at stable (solid lines) and unstable (dashed lines) equilibrium states against phage loss rate constant γ in the Model 3 (case (b)). Bacterial populations $C_{\mathrm{S}}$ (black) and $H_{\mathrm{L}}$ (blue) and the phage population V (red) are shown. Panels (b), (c) and (d) show the critical value of φ necessary to ensure the stability of the lysogenic equilibrium, versus α, η and K. In each case, the parameter space above each curve corresponds to a stable lysogenic equilibrium (LE${}_{\phi }$) while the space under each curve corresponds to stability of the all existing equilibrium (AEE${}_{\phi }$).

4. Discussion

We explore the dynamics of lysogenic phage in a bacterial biofilm, for bacterial hosts both with and without CRISPR immunity. Classical models of phage–bacteria interactions have previously demonstrated that lysogeny can promote the stable co-existence of bacteria and phage [13]. In agreement with these findings, Model 1 explores the baseline conditions under which lysogeny exists in a biofilm, and provides conditions for the stability of a lysogenic equilibrium. The second and third models demonstrate the powerful effect of the CRISPR-Cas system in comparison with the non-CRISPR bacterial population in the first model. The second model predicts that at realistic parameter values, only prophage-free CRISPR-immune bacteria stably survive in a bacterial biofilm. In rare cases, the CRISPR system stably co-exists with the phage, however, the CRISPR system does not allow lysogens to co-exist at this equilibrium. In the third model, there is no CRISPR system present in the lysogens. In the absence of an external source of lysogens (case (a)), CRISPR-immune bacteria are predicted to dominate the biofilm. However if an external source of planktonic lysogens contributes to the biofilm, CRISPR bacteria may co-exist with lysogens or can even be eliminated from the biofilm. This last result is of clinical relevance because CRISPR bacteria are highly resistant to phage therapy. Once CRISPR bacteria are removed from the biofilm, phage therapy has a much higher chance of success.

In particular, case (b) of Model 3 was specifically designed to explore the possibility of using lysogens to penetrate the biofilm and reduce the population of CRISPR-immune bacteria. Unfortunately, the potential success of this strategy is tempered by the results of our numerical work, indicating that large magnitude changes in any one baseline parameter value would be required to eradicate CRISPR-immune bacteria from the biofilm. This suggests that means of varying several parameters simultaneously might hold more therapeutic promise. Although our model treats only one type of virus and corresponding prophage, another interesting possibility is that diverse prophage could be introduced via the lysogens joining the biofilm. In this way, lysogens could produce a number of different viruses via induction, and this could ultimately reduce the entire biofilm population. Since older CRISPR spacers in the bacteria become less efficient, the possibility of CRISPR failure increases with phage diversity, which helps to eradicate CRISPR-immune bacteria $C_{\mathrm{S}}$. Once CRISPR bacteria have been eliminated, the biofilm can be treated by classical therapeutic techniques.

Although the arguments above are highly speculative, the main results of our research are encapsulated in Figure 5 which is divided into four panels (vertical lines). On the left, we use Model 2 to simulate a pathogenic biofilm that is resistant to phage therapy because it consists entirely of CRISPR-immune bacteria. In the second panel, an external, possibly therapeutic source of planktonic lysogens is applied, and the all existing equilibrium state emerges. When the concentration of lysogens in the external source is further increased (3rd panel), the CRISPR-immune bacteria are eradicated from the biofilm and the model shows the same behaviour as Model 1 with only lysogens and phage surviving. At realistic baseline parameters, the lysogenic equilibrium emerges. Once the CRISPR population has been eliminated, the external source of lysogens can be removed and the system is described by Model 1. Variations in a number of parameter values can then be used to eliminate the lysogens (panel 4). For example, as mentioned above, diversity in the prophages carried by lysogens could produce a number of distinct phage populations inside the biofilm, which could infect bacteria and also increase the prophage induction rate. Extending the analysis presented here to include multiple phage types is an intriguing possibility for future work.

Figure 5. Population dynamics of phage-bacteria interaction in Models 2, 3 (case (b)) and 1. Three Vertical lines are drawn to separate the simulations of the models into four panels. Model 2 is simulated for baseline parameters in the first panel, Model 3 (case (b)) is first simulated for baseline parameters and then simulated after an increase in φ in the second and third panels respectively. In the fourth panel, simulations of Model 1 are presented for baseline parameters (top two curves) and for varied parameters based on possible therapy, i.e. $\eta =0.45$, $\alpha =0.3$ and r=0.75, (bottom two curves with increased line-width). The dashed lines represent $C_{\mathrm{S}}$, dash-dot lines represent $C_{\mathrm{L}}$ in Model 2, and $H_{\mathrm{L}}$ in Model 1 and 3, whereas dots symbolize the population density of phage V.

Future work could also relax some of the assumptions and limitations of our model. Most importantly, we were motivated by the possibility of using lysogens to deliver phage therapy to a resistant subpopulation of bacteria, for example within a human host. Biofilms offer an important example of such a subpopulation, but in order to focus on immunity and lysogeny, our model aggregates a spatially complex biofilm population into only two or three compartments, defined by lysogeny and CRISPR status. We thus ignore spatial variations in cell attachment and detachment, cell growth, and spatial variations in the infection rate. A more nuanced spatial model which tracks, for example, the depth of each cell within the biofilm could undoubtedly offer further insight.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC).