Team:Edinburgh OG/Model:Results

PhagED: a molecular toolkit to re-sensitise ESKAPE pathogens

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 phage-bacterial 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) CRISPR-Cas9 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 CRISPR-Cas9 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 two-phage strategy, upon infection of lysogenic bacteria with a lytic phage, the CRISPR-Cas9 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 CRISPR-Cas9 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 CRISPR-Cas9 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 in-text 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·10-8 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 in-text 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 10-15 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 co-exist in chemostat in an equilibrium over a long-time 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 phage-bacteria 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 two-phage re-sensitation 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 10-15 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 re-infection 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 re-infection 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 re-infection 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 re-infection (a2c6858), (C) simulation without superinfection (f192168), (D) simulation without re-infection 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 10-3 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 x-axis 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 10-100 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 r-values 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 r-value	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.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
Half-saturation constant	Km	0.065	1.000	-0.355	-0.995
Half-saturation 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.000	0.634	1.920	1.000
Maximum specific growth rate	μ_I	0.000	1.000	0.000	0.000
Half-saturation 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 r-value. 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 r-value, Ki, μ and b have the highest elasticity values. One per cent increase in Ki will result in 39.75% increase in r-value, meaning the proportion of lytically infected cells will increase dramatically. The main rationale behind the r-value is to quantify the selective pressure that a lytic phage is exerting onto the bacteria. The higher the r-value, 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 r-value (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₀).

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 non-linear. 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 non-linear 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 r-ratio, 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 r-value. 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 pre-selected 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 r-ratio 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).

Model Homepage

Background Methodology

Results References