Non-conventional therapeutic technique to replace CRISPR bacteria from biofilm by inducible lysogen

ABSTRACT Bacteriophage can be an effective means of regulating bacterial populations when conditions allow phage invasion of bacterial colonies. Phage can either infect and lyse a host cell, or insert their DNA into the host cell genome; the latter process is called lysogeny. The clustered regularly interspaced short palindromic repeat (CRISPR) system, linked with CRISPR-associated (Cas) genes, is a regulatory system present in a variety of bacteria which confers immunity against bacteriophage. Studies of the group behaviour of bacteria with CRISPR/Cas systems have provided evidence that CRISPR in lysogenized bacteria can cause an inability to form biofilm. This allows CRISPR-immune bacteria in biofilms to effectively resist phage therapy. Our recent work has described a potential therapeutic technique to eradicate CRISPR-immune bacteria from a biofilm by a continuous influx of lysogens carrying an identical phage sequence. However, this model predicted that the CRISPR-immune population could persist for long times before eradication. Our current focus is on the use of diverse lysogens against CRISPR-capable bacterial populations. The goal of this work is to find a suitable strategy which can eradicate bacteria with a CRISPR system through the influx of finite amounts of distinct lysogens over fixed intervals.


Introduction
Bacteriophage (viruses that infect bacteria) adapt rapidly, such that bacterial defence mechanisms must continually evolve in order to sustain long-lasting phage-bacteria interactions. As each species exerts selective pressure for novel adaptation in the other, a complex co-evolutionary dynamical system emerges. Clustered regulatory interspaced short palindromic repeats (CRISPR) and the CRISPR-associated (Cas) genes constitute one of the advanced defence mechanisms of bacteria. CRISPR/Cas systems not only provide defence against external threats but also regulate the genomic structure of bacteria [7,43,48]. Since CRISPR systems were first discovered a decade ago, their remarkable abilities in genome regulation have been extensively studied, and diverse practical applica- p i Probability of induction tions such as genome editing/sequencing and immuno-therapy have emerged [22,41,49].
In phage-bacteria interactions, CRISPR/Cas systems provide an effective defence against rapidly mutating bacteriophage due to their ability to confer rapidly adapting immunity [6,23]. In the past few decades, increasing concern over antibiotic resistance has revived the idea of using phage therapy, which can potentially eliminate bacterial infections, for medical purposes ( [3,35], but also see [13]). Bacteria containing CRISPR/Cas systems, however, may be immune to a wide range of bacteriophage, and may present a serious obstacle to phage therapy techniques [9,25]. Unfortunately, infections by bacteria equipped with CRISPR/Cas systems are inevitable since estimates indicate that approximately half of bacterial and archaeal species contain these defences [7,12,46]. In addition, CRISPR/Cas mechanisms can provide an even stronger defence to bacteria in biofilms [24], a context in which phage therapies are already very difficult to implement [18,31,37]. Bacteriophage typically infect and kill bacterial host cells to produce more identical phage. However some phage, instead, have the ability to incorporate their genome into the bacterial genome. Bacterial cells containing the viral genome are referred to as lysogens; the viral genome incorporated into the bacterial DNA is called prophage. Lysogens can revert to productive infections at later times, infecting other bacteria. A functioning CRISPR/Cas system, however, does not allow lysogeny in the bacterial genome [14], and in fact very few prophage are found in the genomes of bacteria containing CRISPR/Cas systems [31]. A recent study has provided strong evidence that a particular type of CRISPR/Cas system (type 1F) prevents bacterial cells from forming biofilm, if the cell has been lysogenized by a bacteriophage [21,48]. Another study asserts that the CRISPR/Cas system can induce the death of a lysogen if the lysogen is infected by either the prophage that is induced by temperature increase or a plasmid insertion into the bacteria containing protospacer that is identical to the one present in already existing prophage [14,16,26,45]. In a recent study, CRISPR/Cas system (type III) has shown conditional tolerance over poly-lysogenized bacteria by co-existence, however, the system put their hosts to fitness disadvantage in mixed populations of Staphylococcus aureus [17,38,42].
A number of mathematical models have been developed to study the ecological and evolutionary dynamics of the CRISPR-mediated immune response against rapidly adapting bacteriophage ( [11,19,20,32], for review see [28,29]). The population dynamics of phage, sensitive and resistant non-CRISPR bacteria, and sensitive and immune CRISPR bacteria were first investigated in [32], with a clear focus on experimentally testable predictions, and in the absence of adaptation or mutation. The co-evolution of phage and CRISPR-encoding hosts was incorporated in later models, which predicted the evolution of highly heterogeneous communities, in which newly emerging phage types replace the dominant types by selective sweeps, or previously dominant types reemerge [11]. Of particular interest was the prediction that 'coalitions' of strains with similar resistance phenotypes, but distinct genotypes (bacterial strains with different CRISPR spacers representing the same phage) should be observed. A recent study concludes that the CRISPR/Cas system is responsible for co-evolutionary arms race between bacteria and phages in order to provide survival advantages to both the host and the bacteriophage [44].
Models of evolutionary dynamics in the presence of CRISPR have also demonstrated that selective pressures imposed by phage can explain the empirical observation within a population, newly acquired spacers display more diversity than older spacers [20]. These differences in spacer diversity are also predicted in a stochastic model of bacteria-phage co-evolution [19], which includes the detailed effects of spacer deletion, mutation and recombination. Moreover, bacteria with CRISPR/Cas system can acquire limited amount of spacer in a given amount of time and therefore can be eradicated by high viral mutation rates, suggesting that phage evolutionary behaviour is reciprocal to host immunity adaptability [47].
Building on these previous contributions, we developed an approach to explore the role of CRISPR-encoding bacteria in a biofilm [4]. In particular, three simple models were developed and analysed to predict the impact of CRISPR-immunity on the efficacy of phage therapy against a bacterial biofilm. An important novelty in the approach was to allow lysogens to invade the biofilm, since lysogens without CRISPR systems retain the ability to form biofilm. The lysogens can later be induced to produce viruses to infect CRISPRimmune bacteria. The model was the first approach in which phage-bacteria interactions were modelled within a biofilm that included lysogens, CRISPR-encoding bacteria, or both [4]. This initial model, however, included only a single type of virus, producing a single prophage, whereas the use of multiple viruses is likely to offer a more effective therapeutic technique.
In the sections to follow, we extend this approach to include the possibility that several distinct phage may be used to eradicate a biofilm infection. In the presence of multiple phage, CRISPR bacteria must acquire immunity to every phage to ensure survival. Moreover, in this study, successive as well as simultaneous insertion of distinct lysogens are investigated to develop an optimal therapeutic technique, in which lysogens could successfully invade the biofilm in minimum time.

Model
We consider a model of N distinct phage populations, represented by a vector notation V. The parameters and the time-dependent variables used in this model are defined in (Table 1). Phage adsorb to susceptible cells at rate β, resulting in lysogeny with probability p L or lysis with probability 1 − p L . Lysogeny will increase the density of lysogens carrying the specific prophage, whereas lysis increases the corresponding phage density with burst size b. The phage are lost at rate γ by dissolution or by adsorbing to bacterial cell debris.
Two sets of bacterial populations compose the biofilm, i.e. bacteria carrying prophage L, also known as lysogens, and bacteria with a CRISPR/Cas system C. Each element in these sets represents a bacterial population that grows logistically withthe maximum growth rate r and shared carrying capacity K. Each of these populations leaves the biofilm at rate η, the 'sloughing off' rate. The set of lysogens is indexed by the prophage sequences present in the bacterial genome. Since N is the total number of virus types in the model, there are 2 N types of lysogen. For example, if N = 3 there are eight populations of lysogens, represented by L ijk where i, j and k are either 0 or 1, indicating the absence or presence of a specific prophage in the bacterial genome. Thus, L 000 contains no prophage, whereas L 111 carries all three prophage sequences. Each prophage is induced at the baseline rate α which results in the production of a burst of b bacteriophage and the death of the lysogen. The bacteriophage produced are all of a single type, corresponding the prophage sequence that was induced (see description of matrix M L below).
The second set of bacterial populations C contains CRISPR/Cas systems which act as an adaptive immune system. Each element of this set represents a population that is immune to specific phage types, corresponding to phage sequences that have been previously acquired by the CRISPR/Cas system. Phage whose sequences have not been acquired by a specific cell type can adsorb at rate β. For example, phage j can infect population C 101 . In this case, the cell may acquire the j-type CRISPR-spacer with probability p S , moving to the C 111 population, or lysis may occur. Phage may also adsorb to CRISPR-immune bacteria, that is, there is a small probability of CRISPR failure p F , and in this case lysis of the cell occurs. These assumptions lead to the following general set of equations (the LCV model): where φ(t) is called as a forcing function that represents the input of lysogens L to the system via lysogen therapy. Here, we study the sequential application of a defined sequence of lysogens to the biofilm. Note that each lysogen may contain one or more distinct prophages. Each lysogen in the sequence has an 'interaction time' that is divided into two phases. These phases are separated by an indicator function φ(t), an element of φ(t), as defined in Algorithm section. In the first 'insertion' phase I p , the lysogen is applied to the environment surrounding the biofilm (φ(t) = 0 ) for a specific period of time. In the second 'competence' phase C p , the external source returns to zero (φ(t) = 0) and the lysogens compete with other cells in the biofilm matrix. At the end of the interaction time, the next lysogen in the sequence is applied which follows the same instructions as above. Symbols M VL (t), M VC (t), M L and M V (t) represent matrices associated with the infection of lysogens, the infection of CRISPR bacteria, the production of free virus by induction and the production of free virus by lysis, respectively. The elements of these matrices are presented as the case study in Algorithm section, while the assumptions used for the development of each matrix in the model are described below.
The matrix M VL represents the loss (or gain) in each lysogen population through infection. Thus, the matrix tracks which cell types are susceptible to which viral types, as well as which new lysogens are produced when lysogeny occurs. An additional parameter, p I , allows partial immunity in the model to be tuned. Specifically, prophage sequences with higher indices may confer immunity against viruses with lower indices; while L 000 is infected at rate β by V 1 , V 1 has a reduced infection rate p I β when infecting L 010 . This allows the model to capture co-evolutionary scenarios, for example, in which 'later' viral subtypes may have emerged from, and thus be closely related to, earlier virus populations.
The second matrix M VC similarly represents the loss (or gain), due to infection, of bacterial cells with CRISPR/Cas systems. The entries in this matrix include the spacer acquisition probability p S and CRISPR failure with probability p F . If an infected bacterium does not have the required spacer, it will shift to a new immunity class with probability p S , while if it has the spacer, it is possible that it fails to recognize the protospacer and thus undergoes lysis with probabilityp F .
The matrix M L represents the contribution made by lysogens, through induction, to each viral population. In the case of multi-prophage in a lysogen, the productivity increases while one prophage get production advantage over the others [40]. Probabilities p 1 , p 2 and p 3 etc. allow arbitrary scaling of the induction rate, such that lysogens carrying one, two or three prophages may be induced at different rates; p i gives the probability that a specific prophage is induced, given i total prophage sequences in the genome. Note that if p 1 = p 2 = p 3 , the net result is a linear scaling, in which for example L 111 is induced overall at three times the rate of L 001 . In contrast if p i = 1/i, the induction rate does not depend on the number of prophage sequences (for example if induction occurs due to environmental stress), and upon induction each prophage in the cell has an equal chance of producing the burst of b virions.
Finally, the matrix M V is given by the sum M V = M CV + M LV . The latter matrices represent the contribution to the free phage populations from CRISPR bacteria and lysogens, respectively. The entries in M CV are straightforward, reflecting infections in which a spacer is not acquired, or infections resulting from CRISPR failure. The entries in M LV include any infections that do not result in lysogeny. Once again, the parameter p I is included such that partial immunity against lower indexed viruses may by conferred by higher indexed prophage.

Algorithm
There are variety of ways to define therapeutic procedures; few are defined in Appendix 2, out of which, one procedural type, i.e. constant phases, is described here in detail. Given the total number of distinct prophage sequences (N), numerical solutions to the system of 2 N+1 + N differential equations (Equations (1)-(3)) are obtained using the standard Runge-Kutta integration (ode15s in Matlab, MathWorks Inc.). The parameters and initial conditions are discussed in the following section.
Possible combinations of prophages harboured in the lysogen are called as therapeutic sequences. For N = 3, four therapeutic sequences are made possible, namely, one prophage followed by two simultaneously, denoted by {1, 2}, three prophages each introduced separately, denoted by {1, 1, 1}, two prophages simultaneously followed by one separately, {2, 1}, and finally the introduction of a lysogen containing all three prophage sequences, {3}. The first three sequences are known as successive therapeutic insertions while the last one is called as simultaneous therapeutic insertion of lysogens. The resulted sets of lysogen and CRISPR bacteria can be written in the vector form as The forcing function φ(t) is used in a vector form in modelled Equation (1). Each element of this vector is an indicator function that is represented by φ(t) with same indices as defined for lysogen. Thus, indices describe which lysogen is inserted into the biofilm. Therefore, it can be written as: Note that the number of elements (populations) in the lysogen vector L and CRISPR bacteria vector C is the same. For N = 3, M VL is given by The matrix M VC is the same as M VL , except that the probability of lysogeny is replaced by the probability of gaining a CRISPR spacer, i.e. p Li is replaced by p Si . Thus, have been introduced for simplicity. The matrix M L occurs in the third equation of the LCV model and is the product of two matrices. The first is the binary matrix representing lysogens that carry a specific type of prophage as a 1, and 0 otherwise. The second matrix is a diagonal matrix which provides the probability of inducing specific prophage. The result is: Finally, matrix M V is the sum of two matrices M CV and M LV . The matrix M CV gives the contribution of each CRISPR population to the free virus populations and can be written for N = 3 as Note that the probability of gaining a spacer against the adsorbed phage is associated with those populations that are not yet immune to a specific phage type, whereas the probability of CRISPR-immunity failure is associated with CRISPR-immune populations. Similarly, for N = 3 the matrix M LV can be written as follows:

Parameters and initial conditions
Baseline parameter values have been chosen wherever possible from previous experimental and theoretical studies [2,4,11,15,27,30,33]. The values for the probabilities of CRISPR acquisition (p S = 10 −5 ) and failure (p F = 10 −5 ) are obtained from previous theoretical studies [4,11]. The probability of establishing lysogeny (p L = 0.01) is defined in such a way that one out of hundred cells adsorbing phage genetic material establishes lysogeny [2,30]. As a non-conservative choice for the study of biofilm eradication, we assume that prophage with higher indices confer no immunity to viruses with lower indices, i.e. the probability of infection by a lower indexed virus, p I = 1. Having said that this parameter can be set experimentally by setting phenotype resistant state of latest inserted bacteria [8,34]. The demographic parameters -carrying capacity K, bacterial growth rate r and phage induction rate and loss rates (α and γ ) -are well defined in the literature as discussed in more detail in [4]. However, here we take a burst size of b = 150 phage/cell that is an average burst size of Escherichia coli and Pseudomonas aeruginosa [2]. The baseline value for the adsorption rate is taken as β = 6 × 10 −9 cm 2 / phage / h; sensitivity analysis was performed (see Results section) due to the high degree of uncertainty in this parameter [8]. The parameter value for the sloughing off rate of the biofilm bacteria is chosen as η = 0.01 h 1 [4,15].
The biofilm formation rate from the source of external lysogens applied during therapy is taken to be 0 = 2.78 cells/cm 2 / h. This estimate assumes that a high density of solitary bacteria can be maintained in the medium adjacent to the biofilm. A more detailed analysis of this parameter value is provided in [4,15]. The parameter values are provided in Table 2 while the initial conditions in this section are C 000 = K, all other populations in C = 0, L = 0 and V = 0. Table 2. Baseline values for the parameters used in this study.

Parameter
Value References

Therapeutic sequences
In the simulation results to follow, we take N = 3 and consider a therapeutic protocol in which lysogens containing one, two or three prophage sequences are continuously introduced (input parameter φ(t) in the LCV model) over fixed interval I p = 5. Each insertion phase is followed by a competence phase. We compare therapeutic sequences with the total therapeutic procedure length of 90 h. For example, if a lysogen with single prophage is introduced, then the insertion phase is I p = 5 h and the competence phase C p = 25 h before the next lysogen population is introduced. For the introduction of a lysogen carrying two distinct prophages, the insertion phase remains the same as 5 h whereas the competence phase increases to 55 h. During the insertion phase, the populations of any pre-existing lysogens and susceptible CRISPR bacteria can be drastically reduced in the biofilm. However, we have observed in the previous studies that extinction during the competition phase of either of the two bacteria takes sufficiently long time that 90-h therapeutic procedure is insignificant. Therefore, there is no need to set a bottom threshold value to consider any of the population below that threshold equals to zero. Moreover, increase in the competence phase for few hours has also shown no effect [4].
'Successive' introductions include either one or two prophages per lysogen, whereas 'simultaneous' introductions include all three prophages in the lysogen. Successive introductions are denoted {1,1,1} for the insertion of three distinct prophage one after another, whereas for example {2,1} denotes the application of a single lysogen population carrying two distinct prophage sequences, followed by another lysogen population carrying the third prophage. Figure 1 illustrates typical results predicted for 90-h therapeutic procedure and four possible therapeutic sequences. Note that each new lysogen class (solid line) was preceded by the extinction of already existing lysogen classes. As explained above, the non-conservative assumption has been made that newly inserted lysogen can be infected by already existing phage populations in the biofilm (p I = 1) if their genetic material is not harboured in the inserted lysogen. Conversely, new phage population was also able to start the lytic/lysogenic cycle in the case of their interaction with already existing bacteria.

Results
The distinct phage produced by a newly introduced class of lysogen infected the corresponding class of CRISPR bacteria. This process replaced an already existing class of CRISPR bacteria with a new class that was immune against all existing phage as shown in Figure 1. However, the total population of each class of CRISPR-immune bacteria was less than the preceding class. As a result, after successive therapeutic insertions (Figure 1( a-c)), in a relatively short time, the total population of CRISPR bacteria was reduced far below the carrying capacity K. The introduction of a lysogen with multiple prophages (Figure 1(a-d)) has shown stronger effect than the successive introduction of single prophage (Figure 1(b)). On the other hand, each new class of lysogen outnumbered the previous lysogen classes and quickly approached the steady behaviour as bacteria reached the biofilm carrying capacity. This was desirable because the elimination of previous lysogen classes stops the production of phage types that can no longer infect the bacterial cells with CRISPR-immunity. The colours of the curves are used to distinct each class of lysogen and CRISPR bacteria and are associated as described in the legend: L 000 and C 000 is for grey, L 001 and C 001 is for blue, L 010 and C 010 is for green, L 011 and C 011 is for cyan, L 100 and C 100 is for red, L 101 and C 101 is for pink, L 110 and C 110 is for cream, L 111 and C 111 is for black. Colours for phage populations are blue for V 1 , cyan for V 2 and black for V 3 . Coloured areas under curves of therapeutic lysogens represent the insertion phase whereas light grey areas under curves correspond to the competence phase.
The growth and decay of bacterial classes were divided into two behavioural periods of the system, see Figure 1. During the first period, existing classes of bacteria extinct from the biofilm with the emergence of new classes of CRISPR bacteria and lysogen. The first period began when a new class of lysogen was introduced into the biofilm. These lysogen rapidly increased in the beginning and bounded by the carrying capacity of the biofilm. At the same time, phage population sufficiently increased in order to significantly infect the CRISPR bacteria. Subsequently, a huge amount of lytic infection took place to CRISPR bacteria that increased the phage population in the system. However, due to CRISPRimmune response, a new bacterial class quickly emerged that was immune against the existing phage. It is interesting to note that the emergence of new CRISPR-immune bacteria was due to the acquisition of phage by already existing CRISPR population and not by the bacterial growth. The first period was crucial in the sense that phage population has shown adverse effects on the population of CRISPR bacteria when lysogen was introduced with multiple prophages rather than a single prophage. Thus, CRISPR population significantly decreased in the presence of multiple distinct phage in the biofilm as compared to single phage, see Figure 1(a) and (d).
During the second period, both bacteria competed with each other to grow in the biofilm. The growth of CRISPR-immune bacteria was affected by CRISPR failure against acquired lytic phage adsorption. At the same time, growth of lysogen was affected by the prophage induction. However, the prophage induction indirectly affected the CRISPR bacteria as well by increasing the phage population in the biofilm. As soon as bacteria reached the carrying capacity, their populations became steady. The time duration of this steady behaviour during competition phase did not show any significant effects on the bacterial populations (see e.g. Figure 1(b) and (d)) which concludes that long competition phases does not affect the overall dynamics of the therapeutic procedure.
An important phenomenon has been observed that has shown the significance of early reduction of CRISPR bacteria by introducing lysogen with multiple prophages. The insertion of lysogen with multiple prophages during the first therapeutic phase has shown more reduction in CRISPR bacteria than the insertion of lysogen with multiple prophages during the later insertion phase, see Figure 1(a) and (c). This is because the early introduction of lysogen with multiple prophages (during first therapeutic phase, see Figure 1(c)) increased the second period of lysogen and has shown higher lysogen population in the biofilm when compared with a single prophage in the early inserted lysogen, see Figure 1(a)). Thus, during the second therapeutic phase, the lysogen raised quickly to high value in the first period of growth with all three prophages while shortened the second period of growth. The reason is the non-conservative condition (p I = 1) that has permitted phage adsorption to all kinds of bacteria accept the ones harbouring identical prophage.

Decrease in CRISPR bacteria by optimizing the insertion phase
At the end of each therapeutic insertion phase of lysogen, the total population density of CRISPR-immune bacteria was calculated to determine the optimal procedure time and best lysogen sequence. We examined each of the four possible sequences for 60, 90 and 180 h therapeutic procedures as shown in Figure 2. We observed that for any duration of the insertion phase, the therapies involving the simultaneous insertion of three prophages ({3}) gave the best results, although the combination of a lysogen with two prophages followed by a lysogen with one prophage ({2, 1}) was also effective. Therefore, the results mentioned in Figure 2 can be considered as an extended version of Figure 1. The variation in the duration of therapeutic procedures (60, 90 and 180-h procedures) have not shown any significant difference for their respective sequence.
The time duration for the insertion phase have a strong effect on the eradication of CRISPR bacteria from a biofilm. This can be observed for each therapeutic sequence in Figure 2 when the insertion phase was varied from 0 to 15 h. In this figure, decrease in the CRISPR population can be divided into two stages, depending upon the duration of the insertion phase; the first stage decrease was I p ≤ 7 h while rest of the duration of the insertion phase I p > 7 provide second stage decrease. This is due to the fact that lysogen population started their second behavioural period of growth from time > 7 h, see Figure 1. When the insertion phase duration was greater than 7 h, the lysogen were kept inserted into the biofilm while second behavioural period was started. This helped lysogen to increase their population quickly in the biofilm. However, increasing I p > 15 did not provide further reduction in the CRISPR bacteria, see Figure A6 (column 3), because biofilm reached its carrying capacity. Thus, Figure 2 suggested that I p played a critical role in reducing the CRISPR population. In addition, first phase of lysogen therapy for few hours was enough to replace the existing bacterial class of CRISPR bacteria while the best reduction of CRISPR bacteria was possible when the second period of growth began. This means that if identical lysogen would have been re-inserted during the second growth period, CRISPR bacteria can be significantly reduced in a biofilm as shown in Figure 3. The procedure involving re-insertion of identical lysogen played a critical role in reducing the bacteria to its minimum. In Figure 3, three scenarios were drawn out of which first scenario was same as followed in the above figures. Clearly, it can be observed that early insertion of lysogen for 5 h was not a good strategy to follow. Instead, re-inserting identical lysogen during the second growth period significantly reduced bacteria for all therapeutic sequences. Moreover, the re-insertion of lysogen at the earlier time of second growth period has shown better results than the later time. In comparison between Figures 2 and 3, more than 10 h continuous insertion of lysogen reduced CRISPR bacteria to the same level as 5-h therapy with an optimal re-insertion time.

Sensitivity analyses
Three parameters are analysed in this study to investigate their effects on the reduction of CRISPR bacteria in the biofilm. All other parameters are defined in Table 2. Figure 3 provides an extension to the above result in Figure 2 by comparing the insertion phases I p = 5 and 10 for two parameters, biofilm formation rate 0 and adsorption rate β. The biofilm formation rate, that is, the rate at which solitary bacteria join the biofilm, is a critical dynamical parameter but empirical estimates in the literature are unclear. In this study, we varied the biofilm formation rate, 0 , in the interval [0,3] to see its effects on the eradication of CRISPR-immune bacteria from the biofilm. Figure 4(a) and (c) demonstrates that simultaneous insertion of lysogen ({3}) minimized the population density of CRISPR-immune bacteria. This is consistent with the previous result in Figures 1 and 2. The increase in the insertion phase (from I p = 5 to 10) followed a translation of curve by further decreasing the CRISPR-immune bacteria in the biofilm. The adsorption rate, β, was varied from 10 −10 cm 2 phage −1 h −1 to 10 −7 cm 2 phage −1 h −1 . Among the four therapeutic sequences, simultaneous insertion remained the best strategy. When β is increased from low to high values, the population of CRISPRimmune bacteria showed highly non-linear behaviour around 10 −9 while the population behaviour was smooth before and after that interval. Before this non-linearity, the population remained at the carrying capacity because adsorption rate was so low that it did not significantly affect the CRISPR bacteria. The reason of the non-linearity around β = 10 −9 was again the low adsorption of phage into CRISPR bacteria that was not sufficient for 90h therapeutic procedure. Due to this reason, the first period of growth extended in time while the second period of growth started so late that the procedural time ended before the system reached a steady behaviour. This brought an uncertain state in the therapeutic procedure that was necessary to be addressed to understand the effects of adsorption rate. Despite the slow adsorption rate, the therapeutic lysogen sequences with multi-prophage during their first phase of therapy were least affected by this uncertain output. In addition, at comparatively high adsorption rate, the therapeutic procedures were quite smooth as observed in Figure 1 The increase in the duration of the insertion phase I p partially decreased the nonlinearity interval, however, the non-linearity remains intact for β values around 10 −9 .
Thus, the adsorption rate β was dependent upon the therapeutic procedure time, therapeutic sequence and the duration of I p . For completing the analysis of the parameter set 0 and β, a parametric space was defined between the two parameters for each therapeutic sequence, see Figure 5. The parameter 0 did not show any effect on population of CRISPR bacteria for lower values of β as it was observed in the previous figure. However, it is interesting to observe that a small value of 0 managed to reduce CRISPR bacteria to a reasonably low level while β was high. For β > 10 −9 , significant decrease was observed for all values of 0 , particularly in the case of best strategy {3} when C was dropped down to less than 1 cell/cm 2 . The parameter space could not capture the non-linearity that has been observed in Figure 4 because the parameter space was defined on big step sizes of parameter values. Therefore, the linear behaviour clearly showed that increase in either of two parameters can lead us to a better reduction of CRISPR bacteria. Among the probabilities defined in the model, probability of gaining spacer played a significant role in varying  population of CRISPR bacteria in the biofilm as shown in Figures A1 and 6. Note that least value in the colour bar is less than 0.001 while maximum value is at the carrying capacity. It is well established that the acquisition and interference efficiency of CRISPR/Cas system reduces when phage has high multiplicity of infection [32,39]. Under such assumption, when probability of gaining spacer was increased, the model simulations predicted increase in CRISPR bacteria. For values of p s > 10 −4 , CRISPR bacteria remained at the biofilm carrying capacity for the strategies in Figure 6(a)-(c). However, CRISPR bacteria decrease as the probability of gaining spacer is reduced below 10 −4 . In addition, there is possibility of further reduction in CRISPR bacteria by re-inserting the lysogen during the second growth period as describe in Figure 3.

Discussion
The current study extends a previously studied competition model between lysogens with a single prophage and CRISPR-immune bacteria in a biofilm [4]. The previous model provided a thorough stability analysis of the two competing populations, and demonstrated that high rates of lysogen insertion (biofilm formation rates) were necessary to eliminate CRISPR-immune bacteria from the biofilm. This work has provided a proof of concept in the form of a therapeutic model to show that insertion of lysogens with multiple, distinct prophage can replace CRISPR bacteria from the biofilm in a realistic amount of procedural time and at biologically plausible biofilm formation rates.
The therapeutic model follows a general matrix model for an arbitrary number of distinct prophage sequences, and then present results for three distinct prophages and four possible therapeutic sequences ({1, 2}, {1, 1, 1}, {2, 1} and {3}). Among these sequences, the simultaneous insertion of a lysogen with all three prophages was the most promising therapeutic strategy. These strategies further reduced CRISPR bacteria when lysogen insertion phase was prolonged so that it overlapped the two growth periods. Beside this, two growth periods were defined. It has been observed that the first growth period was only important to introduce lysogen in the biofilm so that virions could reach their peak value. As soon as the virions reached the peak, re-insertion of lysogen became effective. However, the effectiveness held until the biofilm reached its carrying capacity. The qualitative results were robust to sensitivity analyses performed on several parameters and conserved across duration of therapeutic procedures, see Figure 2.

CRISPR bacteria are eliminated faster than lysogens
As a result of the insertion of new class of lysogen into the biofilm, a new kind of bacteriophage is produced and adsorbed by the existing class of CRISPR bacteria. Consequently, a new class of CRISPR-immune bacteria emerges. The previous classes of CRISPR bacteria and lysogens are reduced because the newly emerging phage infection is lethal for both. However, the CRISPR population associated with the previous class decreases more rapidly, as some cells acquire the new spacer and contribute to the newly emerging CRISPR class. In Figure 1, we note that successive insertions {1, 1, 1} and {1, 2} do not initially show this behaviour since the population of CRISPR bacteria is quite high. However, this phenomenon is observed as soon as the population decreases through subsequent insertions. This shows that the effect of successive insertion is slower than a simultaneous insertion.
The relation between the rate of decrease in CRISPR bacteria and lysogens of any class can be analytically investigated. Assume that C i and L i are the ith classes of CRISPRimmune bacteria and lysogen present in the biofilm. A new class of lysogen is inserted into the biofilm L i+1 which produces phage V i+1 that is infectious against the ith class bacterial populations. Thus, the ith class populations monotonically decrease after this time. Considering this as a starting time and taking the time derivative of the ratio between C i and L i , we find d dt A negative outcome of the above relation means that C i decreases more rapidly than L i , though, the concentration of V i is very important since this rate of decrease depends upon the relation α < βp F V i . Thus, higher concentrations of V i will increase the possibility of elimination of bacterial class C i before L i from the biofilm. The above relation is derived for lysogens carrying a single prophage, however, the general relation can be written in a similar way as, α < βp F V i , where V i are the phage produced by lysis ofL i .

Optimal number of prophages per lysogen
The successive insertion of distinct lysogens, each set of lysogen carrying a specific prophage, is less effective than the simultaneous insertion of lysogens (i.e. each lysogen harbouring complete set of prophage), as shown in Figures 1-6. This shows that CRISPR-immune bacteria suppress from the beginning of therapeutic procedure. However, successive insertions suppress CRISPR bacteria in multi-phases depending upon the efficacy of lysogen during each insertion phase, see Figure 1. In Figure 1 Figure 1(b). Thus, simultaneous insertion allows greater diversity of CRISPR-immune bacteria and lesser diversity of lysogen classes as compared to successive insertion. Moreover, the phage density is also reduced in simultaneous insertion at the final procedural time which can significantly reduce the multiple infections to the CRISPR bacteria. Consequently, the events of CRISPR-immune failure decrease [39]. Although the simultaneous insertion shows some evolutionary advantage to CRISPR-immune bacteria, it is not well understood how effective it could be during the therapeutic procedure.
The results above suggested that therapies involving multiple prophages in a single lysogen (simultaneous insertion) were the most promising strategies. However, it was predicted in a previous study that CRISPR bacteria can again dominate in the biofilm if the lysogen therapy stopped before the complete eradication of CRISPR bacteria [4]. Therefore, in the long run, it is very likely that CRISPR bacteria will get immunity against the introduced phages and will dominate the biofilm in case of three prophages. We therefore investigated the long-term effects of single insertions, involving lysogens carrying between 1 and 5 prophages simultaneously, see Figure 7. Informed by the results of Figure 2, the insertion phase was set to 5 h, followed by a long competence phase. The procedural time was chosen to be 3 × 10 6 so that the system shows a steady state for enough long period of time while population density lower than 10 −5 is considered zero. The total population of lysogens, CRISPR bacteria and phage at equilibrium are shown in Figure 7. For lysogens with up to three prophages, the population of CRISPR-immune bacteria re-emerged in the long run, despite the promising results of 90-h therapy as shown in Figure 1. However, the simultaneous insertion of a lysogen population with four or five distinct prophages was able to completely replace CRISPR bacteria.
In this study, the simulations were based on an ODE model with mass action kinetics which assumes a well-mixed system. Such assumptions make difficult to tackle spatial aspects of a physical system and can be misleading. Particularly in case when viruses are produced by burst, not by the process of budding. In the case of viral production by bursting bacteria, hundreds of virions are released at a specific location that is surrounded by the host bacteria (C) which may increase the possibility of co-infection. Besides this, there are other aspects of spatial models that may favour the CRISPR bacteria in the biofilm. For example, once CRISPR bacteria is immune to the existing phage populations, a new class of lysogen needs to overtake the already existing lysogen in order to make possible the interaction between new phages with the CRISPR bacteria. Thus, the above presented results encourage us to make more rigorous analysis of this technique by means of spatio-temporal models in order to foresee the future of this therapeutic technique.

Conclusion
This study investigates a novel therapeutic technique by which 'strong' bacteria (CRISPR bacteria) may be replaced by 'weak' bacteria (lysogens) in a biofilm, in a reasonable therapeutic time. The main point of this study is not to find the best therapeutic strategy, but to show a 'proof of principle' that competition between lysogen cultured in laboratory settings and bacteria living in a biofilm can lead to the design of an efficient therapeutic technique.
Ideally, genome editing could predispose the lysogens for programmed death after a specific period of time or in response to a specific biochemical signal [14,42]. Moreover, the time required for the therapeutic insertion of lysogens and the complete eradication of CRISPR bacteria predicted in this study is reasonable, showing potential for more rigorous analysis of this technique and some initial trials in laboratory settings. Ultimately, lysogen therapy, like phage therapy, may become a promising therapeutic technique for individuals experiencing chronic bacterial infections due to genetic disorders.

Probability factors
Four probabilities are mainly used in this manuscript for model development. Most of them are described in the literature. Here we make an analysis of their behaviour to understand their impact on the modelled system (see Figure A1). All four therapeutic sequences are studied that are possible in case of three distinct prophages in the lysogen. Constant therapeutic procedure type is studied that has 90 h procedural duration. Total population of CRISPR-immune bacteria C is calculated at the end of procedure for each value of underlying probability as shown in Figure A1.
In case of probability of lysogeny, C dramatically decreased from the biofilm when lysogen were introduced into the system, i.e. p L > 0. However, when lysogen contained all kind of prophages in their genome, see therapeutic sequence {3}, reduction in CRISPR-immune bacteria did not depend on the probability of lysogeny. The reason is simply because the lysogen already contained all kind of prophages and therefore did not allow any identical phage infections. The latest insertion of lysogen with a single prophage ({1, 1, 1} and{2, 1}) showed an increase in the CRISPR-immune bacteria for high values of p L . CRISPR bacteria was already immune against existing prophage while lysogen efficiency of inducing a new prophage was reduced due to the establishment of lysogeny from already existing phages. Thus, infections to lysogens from already existing phages reduced the effect of producing effective phage into the biofilm. High values of therapeutic sequence {1, 2} brought CRISPR-immune bacteria to the same level as simultaneous insertion. Two other probabilities of gaining spacer p S and CRISPR-immunity failure p F are also considered in this section. Interestingly, increasing the probability of failure does not show adverse effect on CRISPR bacteria near about the parameter value chosen in this study. On the other hand, probability of gaining spacer is highly sensitive in a sense that it shows a linear increase in the CRISPR population as probability of gaining spacer is increased. This probability is therefore further studied with respect to variation in biofilm formation rate 0 as shown in Figure 5.

Appendix 2
Here I describe four different procedural strategies out of which one therapeutic procedure is explained in the manuscript text. The types of procedures are explained by using Table A1 for the therapeutic sequences based on three number of prophages (np = 3). Figures A3-A6 are made coloured to show the difference among the bacterial classes that are associated to different combinations of bacteriophages. The colour association is shown in Figure A2.

Mixed phases
In mixed phases strategy, each insertion phase (I p ) is followed by competence phase (C p ) while the time for each phase is fixed. Moreover, the phase to phase cycle will be repeated according to the number of prophages present in the latest inserted lysogen. For instance, if I p = 5 and C p = 15 Table A1. Four strategies of therapeutic procedures are defined, namely constant phases, mixed phases, semi-mixed phases and separate phases.
Constant phases L 001 : I P ,C P L 001 : I P ,C P L 011 : I P ,2C P L 111 : I P ,3C P L 110 : I P ,2C P L 010 : I P ,C P L 100 : I P ,C P L 100 : I P ,C P Mixed phases L 001 : I P ,C P L 001 : I P ,C P L 011 : I P ,C P ,I P ,C P L 111 : I P ,C P ,I P , L 110 : I P ,C P ,I P ,C P L 010 : I P ,C P L 100 : I P ,C P C P ,I P ,C P L 100 : I P ,C P Semi-mixed phases L 001 : I P ,C P L 001 : I P ,C P L 011 : 2I P ,2C P L 111 : 3I P ,3C P L 110 : 2I P ,2C P L 010 : I P ,C P L 100 : I P ,C P L 100 : I P ,C P Separate phases L 001 : I P L 001 : I P L 011 : 2I P L 111 : 3I P ,3C P L 110 : 2I P ,3C P L 010 : I P L 100 : I P ,3C P L 100 : I P ,3C P Notes: The therapeutic sequences are shown for three number of prophages in the table. I P and C P stand for the insertion phase and the competence phase. The symbol L ijk represents lysogen with prophage index i, j and k ( = 0 or 1). 0 means the absence of prophage, whereas 1 represents the presence of prophage.

Semi-mixed phases
Semi-mixed phases procedural type is described in detail in the manuscript text. This strategy works similar to 'mixed phases' in case of successive insertions when each lysogen has only one prophage. In case of insertion of lysogen with multiple prophages, the phases are separated from each other. For example, therapeutic sequence {1, 2} with I p = 5 and C p = 15, L 001 would be [I p , C p ] = [5,15] and L 110 would be [I p , C p ] = [10,30]. Population dynamics of two bacteria with CRISPR/Cas system and with lysogen, along three types of bacteriophage population are shown in Figure A4.

Separate phases
In this therapeutic procedure, lysogens are consecutively inserted without any delay of the competence phase. The time of each therapeutic phase depends upon the number of prophages present in the lysogen. If we consider the therapeutic sequence {1, 2} and the same values for I p and C p as in previous cases then L 001 would be [I p , C p ] = [5,0] and L 110 would be [I p , C p ] = [10,45]. Similarly other sequences can be traced out from Table A1.

Constant phases
In this therapeutic procedure, the time for the insertion phase is independent of the number of prophages and, therefore, constant. However, the time for the competence phase depends upon the number of prophages harboured in a given lysogen. Following the previous example of therapeutic phase {1, 2}, L 001 would be [I p , C p ] = [5,15] and L 110 would be [I p , C p ] = [5,30]. Figure A6 presents the population dynamics of bacterial and viral populations present in the biofilm.