Team:TUDelft/Model/Toxin Production


Model of phage infection and toxin production

For PHOCUS to work, administering our engineered bacteriophage should result in the production of sufficient toxins, Cry7Ca1 and short-hairpin RNA (shRNA), to kill the locust, and the production of these toxins should be fast. To determine whether we can achieve these goals, we modelled the temporal evolution of the number of bacteria, phages and produced toxin. The specific aim of this model was to investigate if enough toxin can be produced in the locust gut to kill the locust within 7 days (Q1). Subsequently, we characterised this model using a linear stability analysis and evaluated the sensitivity of our prediction with a parameter sensitivity analysis (Q2). Third, we investigated whether the rate at which bacteria acquire resistance poses a threat against the effectiveness of our biopesticide (Q3).

Model explanation

We adapted the model by Beretta & Kuang that describes bacteriophage infection of one bacterial species (Figure 1) [1].

Figure 1. Graphical representation of phage infection and toxin production model. Bacteria (white cells) replicate over time. In the presence of phages (green), they get infected (turned into infected, red cells). Infected cells produce both toxin and more phages, and lyse. Phages and toxins get degraded over time.

In this model, we distinguish four concentrations: susceptible and infected bacteria, phages, and toxin molecules. We model the development of the number of susceptible bacteria in the absence of phages as logistic growth with growth rate \(\mu\) and carrying capacity C. Susceptible bacteria get infected with rate \(\varepsilon\). In our model, infected bacteria no longer multiply, but they do die with rate \(\lambda\), releasing new phages and toxins as they do so. Phages and toxins are also degraded over time with rates \(\gamma\) and \(\eta\), respectively. Additionally, the system is assumed to be ideally-mixed. The dynamics of this model are described by the differential equations (Eq. 1, Eq. 2, Eq. 3 and Eq. 4):

\begin{equation} \frac{dX_S}{dt} = \mu X_S (1 - \frac{X_t}{C}) - \varepsilon X_S P \tag{1} \end{equation}
\begin{equation} \frac{dX_I}{dt} = \varepsilon X_S P - \lambda X_I \tag{2} \end{equation}
\begin{equation} \frac{dP}{dt} = - \varepsilon X_S P - \gamma P + \beta \lambda X_I \tag{3} \end{equation}
\begin{equation} \frac{dT}{dt} = -\eta T + \alpha \lambda X_I \tag{4} \end{equation}

Model validation

To determine the accuracy of the model, we validated it by applying it to an experimental system consisting of a T7 bacteriophage and E. coli BL21 DE3. To this end, we measured the biomass (OD600) and phage concentration over time. After the cells had grown to the early logarithmic growth phase, phages were added and measurements were made every 5 minutes (see Experiments). The OD600 is a measure of the total bacteria concentration, which includes both susceptible and infected bacteria. We compared the experimental data to the prediction made by our model (Figure 2).

Figure 2. Experimental validation. a) Total cell concentration and b) bacteriophage concentration against time. The red dots are the experimental data and the blue line is the model prediction for the parameter values given in Table 1. We shifted the starting point of the model predictions by 30 minutes to match the experimental data.

Figure 2 shows the results of the measured and simulated concentrations of total cells and phages for phage infection of E.coli using T7 phage. It is to be noted that the experimental concentration of phages (Figure 2b) shows fluctuations, which could be attributed to the fact that in reality a lot of phages adsorb to one cell for a single phage to enter, as only “free phages'' can be measured at any given time. Alternatively, these fluctuations could be explained by the fact that the measurement technique used for determining the phage concentration is prone to human error. The main difference between the model and experimental profiles is that the model underestimates the time delay at which the biomass concentration starts decreasing and the phage concentration starts increasing (0.2h predicted, 0.6h measured). Therefore we introduced a time shift, since t=0 is rather arbitrary, as it is hard to know the induction time in the experiment. The shift suggests that the contact rate constant (ε) might be different from the value we assumed. This is understandable as in reality ε changes with the number of phages i.e. initially ε is very low and needs to cross a “phage replication threshold” before it can become the value we use, which is outside the scope of our model [4]. Nevertheless, the predicted concentration profiles by our model shows clear similarities to the experimental profiles in terms of trend and value. This demonstrates that our equations and parameters are a fair representation of reality.

Q1. Is there enough toxin production to kill the locust within 7 days?

The ultimate goal of our bio-pesticide PHOCUS is to kill the locust within 7 days (Q1). To determine if this goal is feasible, we used the model described above to predict the toxin concentration profile as a function of time.

Figure 3 shows the concentrations of susceptible bacteria, infected bacteria, phage and toxin concentrations of Cry7Ca1 and shRNA over time, as calculated from our simulations.

Figure 3. Phage infection and toxin production prediction in the desert locust gut. a) Susceptible bacteria concentration, b) infected bacteria concentration, c) bacteriophage concentration, d) Cry7Ca1 toxin concentration* and e) shRNA concentration against time. *Where Cry7Ca1 mass is 129kDa [11]. The parameter values for this prediction are found in Table 2.

The simulation shows that within approximately 25 minutes, all susceptible bacteria are infected by the phages (Figure 3a). As the infected cells are lysed, the phage, Cry7Ca1 and shRNA concentrations increase to their maxima (Figure 3c/3d/3e). According to our predictions, PHOCUS will produce a maximum Cry7Ca1 concentration of 1.05 \(\mu\)g/mL in the locust gut (Figure 3d). While it is unlikely that we actually target 100% of the bacteria in the gut, we can use a phage cocktail to ensure that we target sufficient bacterial species. In a study by Wu et al., the Cry7Ca1 concentration required to kill a 50% of a locust population (LC50) within 7 days was already estimated at 0.87 µg/mL [11]. Moreover, the aforementioned LC50 was the concentration Cry7Ca1 administered to the locust and not actually the concentration that reached the gut itself. The advantage of PHOCUS is that the toxin production actually happens locally in the gut and therefore even lower concentrations of Cry7Ca1 could be sufficient to kill the locust.

Similarly, the model predicts a maximum shRNA concentration of \(1.21 \cdot 10^{12}\) shRNA/mL in the locust gut (Figure 3e). Based on some rough calculations, this would result in \(1.57 \cdot 10^9\) to \(6.29 \cdot 10^9\) molecules of shRNA being transferred into the hemolymph (see detailed calculations for more information). A study on RNAi in the desert locust showed that approximately \(4.6 \cdot 10^{10} \) shRNA molecules are required in the hemolymph to induce a significant knockdown effect [10][12] . The two main bottlenecks here are the transport efficiency of shRNA to the hemolymph (1% - 4%, see calculations) and the high degradation of shRNA in the gut. In the design of PHOCUS, we therefore incorporated a strategy to protect and enhance the delivery of shRNA to the hemolymph using special cargo proteins (design page). Moreover, the mode of action of Cry7Ca1 would enhance the delivery of shRNA to the hemolymph as it creates pores in the gut lining of the locust.

One of the design requirements for our biopesticide PHOCUS is that the toxins should be produced fast. The simulation predicts the maximum concentration for both Cry7Ca1 and shRNA is reached after just a few hours. In reality however, the time it takes for the locust to digest the phages should also be taken into account. Additionally, the effective contact rate (\(\varepsilon\)) could be lower in reality as the locust gut does not resemble an ideally-mixed system, which could also slow down the toxin production. Nevertheless, we show that phage infection and toxin production is a rather fast process, which should enable the production of toxin within a day and therefore allow for the locust to be terminated within 7 days.

Q2. How strongly do our predictions depend on the values of the parameters?

In the previous section, we identified some potential bottlenecks of our biopesticide. We therefore aimed to use this model to investigate which design strategies would be most effective in increasing the performance of our biopesticide. To this end, we investigated the behavior of this model using a linear stability analysis (LSA) and investigated the sensitivity of our prediction using a local parameter sensitivity analysis (PSA) (Q2). A LSA can give insight in which parameter changes causes the qualitative behavior of the system to change. The results of this analysis could help us set boundary conditions for some parameters to, for example, force the system into a certain behaviour beneficial to toxin production. Secondly, we performed a local parameter sensitivity analysis (PSA) with the parameters used in Q1 to investigate how changes in these parameters influence our prediction. From this, we could learn which model parameters have the most influence on the toxin production and, thus, which design strategies would be most effective in increasing the performance of our biopesticide.

Linear stability analysis

To better understand the behavior of this system, a common mathematical method called a linear stability analysis (LSA) can be performed. A LSA allows for the characterization of equilibrium points of the system into either linearly stable or linearly unstable. Ultimately, the LSA can give insight into which parameter change(s) cause(s) the stability of an equilibrium point to change, which is known as bifurcation. Knowing the bifurcation points of this system would allow us to better understand the model behavior and set some boundary conditions to certain parameters to, for example, force the system into a certain behavior beneficial to toxin production.

Beretta and Kuang already extensively characterized the equilibrium points of the original model [1]. To adapt the model to our situation, we added a fourth equation (Eq. 4) describing the toxin production (T) profile as a function of time. However, because the original equations describing the susceptible bacteria, infected bacteria and phage concentration are independent of T, the behavior of those equations remains the same.

At an equilibrium point, the time derivatives of all quantities in our model are zero. Our model contains three such points:

\begin{equation} E_{-} = (0,0,0,0) \end{equation}
\begin{equation} E_{+} = (C, 0, 0, 0) \end{equation}
\begin{equation} E_{*} = (\frac{\gamma}{\varepsilon (\beta-1)}, \frac{\gamma \mu (C \beta \varepsilon - C \varepsilon - \gamma)}{\varepsilon(\beta - 1)(C \beta \varepsilon \lambda - C \varepsilon \lambda + \gamma \mu)}, \\ \frac{\lambda \mu (C \beta \varepsilon - C \varepsilon - \gamma)}{\varepsilon (C \beta \varepsilon \lambda - C \varepsilon \lambda + \gamma \mu)}, \frac{\alpha\gamma\lambda \mu (C \beta \varepsilon - C \varepsilon - \gamma)}{\varepsilon\eta(\beta - 1)(C \beta \varepsilon \lambda - C \varepsilon \lambda + \gamma \mu)} \end{equation}

We can determine the linear stability of each equilibrium point by calculating the eigenvalues of the linearized system near that point. If any of the eigenvalues has a positive real part, the solution will tend to move away from the equilibrium point, making the point unstable.

One of the eigenvalues of the linearized system around \(E_{-}\) is equal to \(\mu\), which is a positive number, so \(E_{-}\) is an unstable point. Intuitively this makes sense, as one can imagine that when there is a small concentration of susceptible bacteria, these are able to multiply at an exponential rate thus moving away from this point.

We found that the equilibrium point \(E_{+}\) is linearly stable when

\begin{equation} 0 <\beta < 1 + \frac{\gamma}{C\varepsilon} \end{equation}

However, when

\begin{equation} \beta > 1 + \frac{\gamma}{C\varepsilon} \end{equation}

\(E_{+}\) is unstable and \(E_{*}\) was found to be feasible under these conditions. The linear stability analysis shows that for increasing values of \(\beta\) the solution near \(E_{*}\) bifurcates from a bacteria-phage equilibrium toward a periodic solution.

Thus, this analysis shows that, depending on the viral replication factor (\(\beta\)) of our engineered phage, our system can exhibit three qualitatively different kinds of behaviors: susceptible bacteria outgrow the phages, bacteria and phages reach an equilibrium state and phages kill all susceptible bacteria. For our biopesticide, the last qualitative behavior where all susceptible bacteria are killed is most interesting to produce the highest toxin concentration.

Parameter sensitivity analysis

To date, no research has been done to determine the parameter values used in the model of host-phage infection and toxin production in desert locust gut conditions. Therefore, estimates of these parameter values were derived from literature. Due to uncertainty in these estimates, it was assessed whether the model is locally sensitive to changes in parameter values. An important measure for the effectiveness of our biopesticide is the maximum amount of produced toxin, i.e. maximum Cry7Ca1 and shRNA concentration. Therefore, we performed two local parameter sensitivity analyses (PSA), one with the parameters of the Cry7Ca1 toxin and one with the parameters for shRNA with their respective maximum concentration as output variables. A second important measure for our biopesticide is the time it takes to reach the maximum concentration, as we need our biopesticide to be fast. Therefore, we also assessed the local sensitivity of both the Cry7Ca1 and shRNA parameter sets to the time it takes to reach the maximum concentration. During the sensitivity analysis we varied the parameters between 20% and 500%. We chose this range as some of the parameters, such as the viral replication factor, are known to vary within this range [9][13]. Furthermore, degradation rates, such as phage and toxin degradation, are difficult to predict and therefore a rather broad range was chosen for the sensitivity analysis.

The PSA helped us understand how each parameter affects the system and gives us insight on what could be improved. We showed that the toxin is produced very quickly in most cases and, in terms of speed, no improvements need to be made. The PSA showed us that the amount of toxin produced is sensitive to particular parameters, which we can use to formulate design strategies that ensure that the maximum possible toxin is produced and reaches its site of action.

We learned that increasing the number of target bacteria, increasing the number of toxins produced per cell and increasing the stability of shRNA in the locust gut are the best strategies to improve the effectiveness of PHOCUS. We therefore propose the following design strategies to improve the toxin production:

  • Targeting more bacteria enhances the maximum amount of toxin produced. The best way to target more bacteria is to use a cocktail of phages to attack a majority of species in the locust gut.
  • The Pinetree simulations from UT Austin suggest that rather than replacement of non-essential genes, the insertion of our genes close to well-expressed essential genes, could lead to a 26 fold increase in expression in E.coli with the T7 phage (see Partnership). This result demonstrates the potential of phages to produce high toxin quantities.
  • Strategies to protect the shRNA from degradation would greatly increase the delivery of shRNA to the site of action. Therefore, we designed a strategy that uses a RNA binding domain to protect the shRNA from degradation by RNAses [15]. Increasing the stability of Cry7Ca1 would also be a valid strategy to improve the toxin production. However, we suggest to first focus on improving the stability of shRNA as this would be the most effective way to improve our biopesticide.
  • Strategies to enhance the transport of shRNA to the hemolymph would also greatly increase the delivery of shRNA to the site of action. This is already being done to great extent by the Cry7Ca1 as it’s mode of action is to form a pore in the locust gut membrane. We also have a secondary approach to assist transport using the trypsin modulating oostatic factor [10]. While the modelling of these strategies is out of the scope of this model, they are positive prospects for increasing the delivery of shRNA to the hemolymph of the locust.

Further information of our strategies can be found in the design page.

Q3. What is the effect of bacterial resistance on toxin production?

As a defense mechanism against phage attack, bacteria acquire resistance against phages. Mechanisms include spontaneous mutations, restriction modification systems, and adaptive immunity via the CRISPR-Cas system [16]. PHOCUS depends on bacteria that are susceptible to phage infection to produce our toxin molecules. If the target bacteria is resistant to the phage, the toxins are not produced and the desert locusts will not die. Therefore, resistant bacteria may pose a threat to our project.

We built upon our previous findings by adding an additional layer of complexity to the model. Formerly, we assumed the bacteria population in the desert locust gut to consist of susceptible and infected bacteria, while in reality, bacteria acquire resistance against PHOCUS (Figure 8) .

Figure 8. Graphical representation of phage infection and toxin production model including phage resistance. Suceptible bacteria (white cells) and resistant bacteria (blue cells) replicate over time. In the presence of phages (green), susceptible cells get infected (turned into infected, red cells). Infected cells produce both toxin and more phages, and lyse. Phages and toxins get degraded over time. In time, susceptible bacteria also turin into resistant bacteria.

To investigate the effect of bacteriophage resistance on the maximum toxin production and toxin profile, we described the change in concentration of resistance bacteria over time as logistic growth with growth rate \(\mu\) (equal to susceptible bacteria) and a rate, \(\phi\), at which susceptible bacteria become resistant.

\begin{equation} \frac{dX_r}{dt} = \mu X_r(1- X_t/C) + \phi X_s \tag{5} \end{equation}

Figure 9 shows the concentrations of susceptible, infected and resistant bacteria, phage and toxin concentrations of Cry7Ca1 and shRNA over time, with the initial concentration of resistant bacteria set at 0.

Figure 9. Simulation of the effect of bacterial resistance of phage infection and toxin production assuming no initial resistant bacteria. a) Susceptible bacteria concentration, b) infected bacteria concentration, c) resistant bacteria concentration, d) bacteriophage concentration, e) Cry7Ca1 toxin concentration* and f) shRNA concentration against time. The concentrations were simulated for 10 hours using the parameters listed in Table 3. *Where mass is 129kDa [11].

Here, we see that, in the timeframe that toxins are produced, the concentration of resistant bacteria is negligible (Figure 9c). Thus, even though the LSA indicates that resistant bacteria will take over the gut eventually, the rate at which bacteria obtain resistance is too slow to negatively affect the maximum toxin concentration and toxin profile.

Continue to our 2D biofilm simulation!