Model: results
Assumptions
The model presented below is based on the following assumptions:
 (1) Each bacterium, regardless of the phenotype, can be infected with only one phage at a time.
 (2) Both lytic and lysogenic bacteria can be either reinfected with the same phage or superinfected by any other phage.
 (3) Reinfection or superinfection has no impact on the bacterial life cycle.
 (4) The latency time is constant for every phagebacterial interaction.
 (5) The burst size is constant for each lysis for both phages.
 (6) The specific growth rate of the bacteria is described with the Monod equation, which also depends on the concentration of the growth limiting substrate.
 (7) The induction rate of lysogenic bacteria is constant.
 (8) CRISPRCas9 system encoded by a temperate phage is 100% efficient and upon infection irreversibly renders susceptible cells immune to lytic infection without a time delay.
Model of a continuous process
Here we present a new model describing interactions of one species of bacteria and two phages: temperate and lytic. The list of interactions is described as follows: the initial population of susceptible bacteria (X_{S}) are infected with a population of temperate phage (PT) at time t=0. Bacteria that are get infected with this phage form a new population of lysogenic bacteria (X_{L}). These bacteria obtain protection from a lytic phage (P) and can pass it to their daughter cells upon division. Upon induction at rate q lysogenic bacteria enter the lytic cycle, which leads to the production of viral particles and cell burst after the time latency (T) has passed. After some time (ε), the lytic phage is added into the system. Lytic phage can also infect susceptible bacteria and form a new population of infected bacteria (X_{I}). After the time latency has passed, infected bacteria burst open releasing b phage particles that can then repeat the cycle. Lytic phage can also infect a population of lysogenic bacteria, however upon entering the cell, phage genome gets destroyed by the CRISPRCas9 system delivered by the temperate phage. Similarly, a temperate phage can infect a lytically infected bacteria before it gets burst, however this has no effect on the bacteria as will be discussed further (Figure 1).
Figure 1. Mapping of interactions between phages and bacteria in a new model. X_{S} – susceptible bacteria, X_{L} – lysogenic bacteria, X_{I} – bacteria infected with a lytic phage, P – lytic phage, PT – temperate phage.
Models by Levin et al. (1977) and Qiu (2007) set the foundation for a new model presented below. The consumption of resources is described as before in both models using Monod equations by all populations of bacteria.
Population of susceptible bacteria decreases as it gets infected by either of the phages proportionally to their absorption rate constants (Ki for a lytic phage and KiT for a temperate phage) and the concentrations of phages and bacteria.
Population of lytically infected bacteria increases upon infection of susceptible bacteria with a lytic phage and decreases as cells infected T time ago lyse.
The population of lysogenic bacteria rise with the infection of X_{S} with a temperate phage or via normal cell division similarly to X_{S}. As described in Qiu (2007) these bacteria can be induced at a rate q to leave the lysogenic cycle and enter the lytic cycle leading to cell lysis. We also incorporated the removal of induced bacteria from the chemostat before lysis as well as the time latency (T_{T}) between induction and lysis using DDEs as described in Levin et al. (2007).
The population of lytic phages increases as lytically infected bacterial cells T time ago lyse, releasing the phages. The population of lytic phage can infect all populations of bacteria, of which only X_{S} will result in an outcome by rendering cells into X_{I}. To highlight the essence of the proposed twophage strategy, upon infection of lysogenic bacteria with a lytic phage, the CRISPRCas9 system inside lysogenic bacteria cuts the lytic phage genome, thereby protecting the bacteria from its viral genes.
The population of lysogenic phage increases proportionally to the induction rate and the concentration of lysogenic bacteria that were induced T_{T} time ago. Similarly to P, PT is able to infect all
populations of bacteria, from which only X_{S} can be turned into X_{L}. Although all lytic phage genomes carry protospacers that can be cut by the CRISPRCas9 system carried by the temperate phage,
we argue that upon superinfection of X_{I} with PT, in most cases there will not be enough time for the temperate phage genes, including CRISPRCas9 system, to be translated and produce a significant enough
number of protein molecules to completely stop or divert the lytic cycle to a lysogenic one.
b b_{T} 
98.0 (Levin, Stewart and Chao, 1977) 
P 
1.0·106 particles (Levin, Stewart and Chao, 1977) 
C 
7.0·109 cells (see intext calculations) 
q 
0.35 per hour (Qiu, 2007) 
D 
0.2 per hour (Levin, Stewart and Chao, 1977) 
S 
30.0 μg/ml (Levin, Stewart and Chao, 1977) 
d 
0.1 per hour (Suttle and Chen, 1992) 
T T_{T} 
0.5 hours (Levin, Stewart and Chao, 1977) 
Ki K_{iT} 
6.24·108 ml/h (Levin, Stewart and Chao, 1977) 
X_{S} 
1.0·104 cells/ml (Levin, Stewart and Chao, 1977) 
K_{m} 
4 μg/ml (Levin, Stewart and Chao, 1977) 
X_{I} 
0 
μ_{max} 
0.738 per hour (Levin, Stewart and Chao, 1977) 
Y 
3.85·105 (Levin, Stewart and Chao, 1977) (see intext calculations) 
Table 1. Initial variable values for the simulations.
Figure 2 shows a simulation of a new model outlined above with the parameters listed in Table 1. Simulation starts at t=0 with the initial X_{S}= 1.0·104 cells/ml and PT= 1.0·106 particles/ml, and both X_{I} and X_{L} set to 0. As time progresses, susceptible bacteria get infected with PT forming a new population of X_{L}. At t=5 hours a lytic phage is introduced into the system at a starting concentration of 1.0·104 particles/ml and starts to infect X_{S}, converting them to X_{I}. At t=5.5 hours the first population of X_{I} lyses, with each bacterium releasing 98 phage particles into the solution, which repeat the cycle. At about t=10 hours, the rate of infection of X_{S} by populations of both phages, combined with the rate of outflow of cells due to dilution, becomes equal to the growth rate of X_{S}, leading to the levelling off of the solid red line representing X_{S}. This is followed by a further slowly increasing decline of X_{S} due to the increase in population of both phages. At about t=20 hours the concentration of X_{S} drops below 1015 particles/ml, leading to extinction. Meanwhile, the concentration of X_{L} keeps growing steadily fuelled by the conversion of X_{S} to X_{L} up until about t=13 hours. It then continues a slower growth due to depletion of X_{S}, but fuelled by a steady natural cell division. Figure 3A shows simulation of the same model over 250 hours. It shows that in this particular simulation populations of temperate phage and lysogenic bacteria can coexist in chemostat in an equilibrium over a longtime period.
Figure 2. Simulation of the new model in a chemostat. All parameters set as in Table 1, ε = 5 hours (d13160f). X_{S} – concentration of susceptible bacteria, X_{L} – concentration of lysogenic bacteria, X_{I} – concentration of bacteria infected with lytic phage, P – concentration of lytic phage, PT – concentration of temperate phage, S – concentration of substrate.
For comparison, we simulated all models that we used for the development of a new model using the same set of initial parameters (Table 1) next to each other (Figure 3). The plot of X_{S} is somewhat similar in Figure 3B and C, whereas the dynamics of X_{S} in Figure 3A and D is completely different from the two before and from each other.
Figure 3. Simulations of the models of phagebacteria interactions in a chemostat. (A) New model (1f6ea35), (B) model by Levin et al. (1977) (5257aa0), (C) model by Campbell (1961) (c8eebe6), (D) model by Qiu et al. (2007) (f001cb2).
Model of a batch process
The model of a batch process for our set up was developed from the model of a continuous process described above by removing the inflow of nutrients and outflow of culture from all equations. The final system of DDEs is as follows:
One of the key questions for the twophage resensitation approach is to define the time when the population of antibiotic resistant bacteria, here referred as susceptible bacteria, goes extinct. Therefore, in all simulations we drew a particular attention to the time at which the concentration of X_{S} goes below 1015 particles/ml as explained in Methodology. In addition, we calculated the ratio of average concentration of lytically infected bacteria to the average concentration of susceptible bacteria (r) also as described in Methodology.
In order to see how inclusion of the reinfection and superinfection affects the time of extinction of X_{S}, in Figure 4 we simulated a full model and variations with and without these two processes. It is clear that with the inclusion of the superinfection has little or no impact on the extinction time (Figure 4A and C). Whereas exclusion of the reinfection process from the model results in decreased extinction time. Although from these simulations it might be suggested to exclude the superinfection from the final model for the further simulations, we decided to leave due to the following reasons: (1) the absence of a noticeable effect of superinfection might be just a result of the simulation with a given set of initial parameters and might change together with the parameters. (2) we could not find any previous studies of phage and bacteria models with either superinfection or reinfection included into the process, so it is an opportunity to present a new data here.
Figure 4. Simulations of final batch models with parameters mentioned in Table 1. (A) A new complete batch model (1c62895), (B) simulation without reinfection (a2c6858), (C) simulation without superinfection (f192168), (D) simulation without reinfection and simulation (84e0f30). The time above each plot shows the extinction time of X_{S}. X_{S} – concentration of susceptible bacteria, X_{L} – concentration of lysogenic bacteria, X_{I} – concentration of bacteria infected with lytic phage, P – concentration of lytic phage, PT – concentration of temperate phage, S – concentration of substrate.
Another important variable of the batch model is the time at which the lytic phage is added to the model. Although the variation of this parameter in the continuous model showed no impact on the extinction time, we expected to see the difference in the batch model. Interestingly, the batch model also showed little variation in the output up if the phage is added at any time up until 9 hours after the start of the simulation (Figure 5). Figure 5D shows that introduction of the phage at 9 hours sharp results in the lowest extinction time (10 hours). However, from the same plot it is also seen that the population of lysogenic bacteria levels off at the concentration around 103 particles/ml, which should not happen due to constant lysis of these cells. This issue is further investigated in sensitivity analysis.
Figure 5. Simulations of a batch model with different times of addition of lytic phage. (A) P added at t= 0 (91b1128), (B) t=5 (8fbc2da), (C) t=7 (6376bb5), (D) t=9 (eaf2859), (E) t=11 (863e5aa), (F) without P (1aa3048). Red arrow below xaxis shows when P is added to the model, r – ratio of average concentrations of X_{I} to X_{S}, t – time when X_{S} goes extinct.
Sensitivity analysis
All previous simulations of the continuous and batch models assumed the parameters, such as the burst size, phage absorption rate and time latency, etc. to be the same for both lytic and temperate phage as listed in Table 1. Similarly, there was no difference in the parameters of all three populations of bacteria: susceptible, lytically infected and lysogenic. Biologically, it is unlikely that all of these parameters will be the same across different populations and species. Therefore, the initial model of the batch process was extended to take allow variation in all 20 parameters.
We then performed a local sensitivity analysis using the expanded model, each parameter was changed from its nominal value (Table 1) by either subtracting 50% of the nominal value to obtain the minimal value or by adding 50% to obtain the maximal value. The range between two values was then divided into 10100 equal parts and the simulation was then run each time taking a new value from the range for the parameter in interest. Both the extinction time of X_{S} and rvalues were tracked in every run. This data was then used to calculate the elasticity of the parameter with respect to the two outputs and the Pearson Correlation Coefficient (PCC) as described in Methodology.
Parameter 
E of X_{S} extinction time 
PCC 
E of rvalue 
PCC 

Maximum specific growth rate 
μ 
0.398 
0.995 
28.800 
0.976 
Maximum specific growth rate 
μL 
0.231 
1.000 
0.977 
0.946 
Burst size 
b_{T} 
0.190 
0.982 
0.978 
0.851 
Starting concentration of substrate 
S_{0} 
0.130 
0.951 
0.186 
0.935 
Starting concentration of susceptible bacteria 
X_{S0} 
0.129 
0.990 
0.350 
0.998 
Adsorption rate 
Ki_{T} 
0.112 
0.991 
0.997 
0.737 
Induction rate 
q 
0.110 
0.943 
0.841 
0.892 
Halfsaturation constant 
Km 
0.065 
1.000 
0.355 
0.995 
Halfsaturation constant 
Km_{L} 
0.037 
0.999 
0.439 
0.997 
Starting concentration of temperate phages 
P_{T0} 
0.019 
0.999 
0.948 
0.924 
Bacterial yield 
Y 
0.004 
0.949 
0.000 
0.000 
Bacterial yield 
Y_{L} 
0.003 
0.775 
0.000 
0.000 
Latency time 
T_{T} 
0.003 
0.723 
1.265 
0.988 
Burst size 
b 
0.001 
0.792 
15.429 
0.992 
Adsorption rate 
Ki 
0.000 
0.066 
39.750 
0.976 
Latency time 
T 
0.000 
0.441 
0.600 
0.202 
Starting concentration of lytic phages 
P_{0} 
0.000 
0.634 
1.920 
1.000 
Maximum specific growth rate 
μ_{I} 
0.000 
1.000 
0.000 
0.000 
Halfsaturation constant 
Km_{I} 
0.000 
0.999 
0.000 
0.000 
Bacterial yield 
Y_{I} 
0.000 
0.960 
0.000 
0.000 
Table 2. Results of parameter sensitivity analysis on the time at which X_{S} goes extinct. Elasticity (E), Pearson Correlation Coefficient (PCC). Parameters of susceptible bacteria and lytic phage are without subscript, L – lysogenic bacteria, I – lytically infected bacteria, T – temperate phage (ae9f326).
The results of the sensitivity analysis are presented in Table 1 and sorted in descending order of the absolute elasticity values of the X_{S} extinction times. To better represent the relative elasticity values of different parameters, they are also presented in spider diagrams (Figure 6). These results clearly highlight parameters that have a significant impact on the outputs: the extinction time of X_{S} and rvalue. The top 3 parameters that have the highest influence on the extinction times are μ, μL and b_{T}. Theн all have a negative linear correlation with the output, for example one per cent increase in μ will result in 0.398 per cent decrease in the extinction time. The relationship is linear since the absolute value of PCC is very close to 1.0.
In terms of the rvalue, Ki, μ and b have the highest elasticity values. One per cent increase in Ki will result in 39.75% increase in rvalue, meaning the proportion of lytically infected cells
will increase dramatically. The main rationale behind the rvalue is to quantify the selective pressure that a lytic phage is exerting onto the bacteria. The higher the rvalue, the larger the proportion of susceptible bacteria that
gets infected with the lytic phage and dies due to lysis. Therefore, one has to be careful whenever trying to manipulate the output variable, such as the extinction time of XS, as an undesirable side effect might result in a significant increase in
the selective pressure thereby raising the risk of development of the resistance to the phage in bacteria.
Figure 6. Spider diagram of the sensitivity analysis results. (A) The effect of parameter change on the X_{S} extinction time (c342063), (B) the effect of parameter change on rvalue (1be38af). Only parameters with absolute elasticities over 0.01 are shown.
Based on the data from Table 2 we selected five parameters that either had significant elasticity of the extinction time or a low PCC to study their change and relationships with the output values in more details: the burst size of both phages (b, bT), induction rate (q), time latency of a temperate phage (T_{T}) and the initial concentration of the lytic phage (P_{0}).
The results of the sensitivity analysis of separate parameters are simulated in Figure 7. It shows that the relationship between the induction rate and the extinction time is nonlinear. This result was not obvious from the PCC value in Table 2, since the percentage change of the nominal value of 0.35 only analysed the lower end of the curve displayed in Figure 7A. The relationship between the time latency of a temperate phage and the extinction time is also nonlinear and shows an obvious minimum extremum of around 0.5 hours at which the extinction time takes the lowest value (Figure 7B). The burst size of a temperate phage is shows a more linear relationship with the extinction time compared to the burst size of a lytic phage, yet the latter has a much larger elasticity of the rratio, confirming its importance in spreading the infection to other cells (Figure 7C and D). There is nearly a linear relationship between the initial concentration of the lytic particles and the extinction time and rvalue. It is worthy noting that this parameter can be easily manipulated by the experimenter, so its influence on the selective pressure has to be always taken into account.
In addition to the five preselected parameters, we also simulated the the variability in the time when the lytic phage is added into the system (Figure 7F). The plot shows that the extinction time is independent from the time of introduction up until 9 hours after the start of the simulation. Although this result matches the one we obtained previously (Figure 5), the pattern of the plot and unexpectedly round number suggest a numerical error during the simulation. So far we were unable to identify the source of this issue, however we suspect it originates from the PyDDE algorithm for solving DDEs.
Figure 7. Sensitivity and relationship analysis of the separate parameters. The extinction time of X_{S} and rratio are plotted against (A) induction rate (d4f3bd5), (B) time latency of a temperate phage (7a2891f), (C) burst size of a temperate phage (c83187d), (D) burst size of a lytic phage (32e0d55), (E) initial concentration of a lytic phage (particles/ml) (df83949), (F) time at which lytic phage is introduced into the system (hours) (eac0ccf). Horizontal line shows the nominal parameters.
Development of stochastic model
We spent a lot of time trying to convert our existing deterministic model into a stochastic one. Most of these attempts failed at the attempt to interpret or replicate a limited number of previous papers that require a much better knowledge of mathematics, which we lacked as a team (Carletti, 2002; Smith and Trevino, 2009). Those few successful attempts produced more or less sensible results, but we still struggled to understand the rationale behind and thus decided to not include them here. Therefore, we made a collaboration with the team from Newcastle. They helped us to develop a basic stochastic model based on Campbell (1961) (Figure 8).
Figure 8. Stochastic simulation of Campbell’s (1961) model. Copy number reflects the total number of particles in the solution, as opposed to concentration in deterministic model. Code is developed in collaboration with Newcastle iGEM team 2017 (1de45de).