Team:TUDelft/Model/Toxin Production

PHOCUS

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}
Details of the model
Susceptible bacteria concentration

The rate of change in the number of susceptible bacteria (\(X_s\)) is the difference between logistic growth (first term) and decrease due to phage infection (second term). Growth depends on the growth rate (\(\mu\)) and the susceptible bacterium concentration and is increasingly limited as the total bacterium population size (\(X_t = X_s + X_i\)) approaches the carrying capacity (C). The carrying capacity represents the maximum number of cells the system can sustain. The rate at which infection of susceptible bacteria occurs depends on the susceptible bacteria concentration, the phage concentration (\(P\)) and the effective per bacteria contact rate constant with viruses (\(\varepsilon\)).

\begin{equation} \frac{dX_S}{dt} = \mu X_S (1 - \frac{X_t}{C}) - \varepsilon X_S P \tag{1} \end{equation}
Infected bacteria concentration

The rate of change in the number of infection bacteria (\(X_I\)) over time is assumed to only increase as susceptible bacteria get infected (first term) and decrease as the infected bacteria are lysed according to a lysis rate (\(\lambda\)) (second term).

\begin{equation} \frac{dX_I}{dt} = \varepsilon X_S P - \lambda X_I \tag{2} \end{equation}
Phage concentration

The rate of change in the number of phages (\(P\)) over time decreases when either phages enter their host bacteria to replicate (first term) or due to phage degradation (second term), and increases when infected bacteria lyse and new phages are released (third term). Phage degradation depends on the phage concentration and a fixed phage degradation rate (\(\gamma\)). The release of new phages is described by multiplying the rate of lysed host bacteria by a constant called the viral replication factor (\(\beta\)).

\begin{equation} \frac{dP}{dt} = - \varepsilon X_S P - \gamma P + \beta \lambda X_I \tag{3} \end{equation}
Toxin concentration

We assumed that the toxin concentration only increases when a cell is lysed and toxin is released into the gut (second term) and decreases as toxin is degraded (first term). Toxin production occurs through the same mechanism as phage production: toxins are produced by the host bacterium due to infection by a genetically engineered phage and toxins are released when the bacterium bursts. Hence, toxin production is described in a similar way as phage production: by multiplying the concentration of lysed host bacteria by a constant that denotes the amount of toxin produced per cell (\(\alpha\)). Toxin degradation depends on the toxin concentration and the overall degradation constant (\(\eta\)).

\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.

Details of experimental validation

We solved the ODEs of the model numerically using the SciPy ordinary differential equation solver for Python (see our Python code). The initial conditions required for solving the system of ODEs were set to the initial conditions of the wet lab experiment. A complete list of the used parameters including sources is given in Table 1.

Table 1. Parameters and initial conditions for experimental validation.
Parameter Value Explanation Source
\(\mu\) 3.0 (1/h) Growth rate of E. coli BL21 DE3 [2]
C \(5.6 \cdot 10^9\) (cells/mL) Carrying capacity (maximal bacteria population) [2]
\(\varepsilon\) \(2.79 \cdot 10^{-9}\) (mL/h) Effective per bacteria contact rate constant with viruses [1]
\(\lambda\) 3.57 (1/h) Bacterial lysis rate [3]
\(\beta\) 179 (-) Viral replication factor [3]
\(\gamma\) 20 (1/h) Phage degradation rate Based on the experimental results
\(X_{s,0}\) \(3.68 \cdot 10^{8}\) (number/mL) Initial concentration of susceptible E.coli cells Based on the experiment
\(X_{i,0}\) 0 (number/mL) Initial concentration of infected E.coli cells Based on the experiment
\(P_0\) \(3.0 \cdot 10^{6}\) (number/mL) Initial phage concentration Based on the experiment

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.

Details of model prediction

We solved the ODEs of the model numerically using the SciPy ordinary differential equation solver for Python (see our Python code). We assumed that the microbiome of the locust gut is in equilibrium and the bacteria have grown to their carrying capacity (C) , so the initial concentration of susceptible cells is equal to the C. For C we assumed the known cell concentration of the locust gut [5]. As initial phage concentration we used \(3 \cdot 10^5\) number/mL. Furthermore, since Cry7Ca1 will be produced in the gut of the locust, we assumed the degradation rate to be equal to that of another Cry toxin at the temperature and pH of the gut [6]. Unfortunately, we could not find a better estimate on the rate of degradation and the value we found does account for the proteolytic activity of the gut. However, Cry toxins have evolved to work in the gut and, therefore, there are some natural mechanisms that protect them [7]. Lastly, for some of the phage and bacteria parameters (\(\lambda\), \(\beta\), \(\mu\)) we generalised these parameters according to literature to match values for average bacteria and phages in gut like conditions. A complete list of parameters and initial conditions is given below (Table 2).

Table 2. Parameters and initial conditions for model prediction.
Parameter Value Explanation Comment Source
\(\mu\) 0.172 (1/h) Bacterial growth rate The growth rate in the gut is much lower, due to anaerobic conditions. Average bacterial doubling time in the mouse gut is assumed to be a good estimate. [8]
C \(3.34 \cdot 10^9\) (cells/mL) Carrying capacity (maximal bacteria population) - [5]
\(\varepsilon\) \(2.79 \cdot 10^{-9}\) (mL/h) Effective per bacteria contact rate constant with viruses - [1]
\(\lambda\) 0.952 (1/h) Bacterial lysis rate λ(SS68C) is a good average for lysis rate as this lysis rate is the slowest lysis rate for a fast phage. [9]
\(\beta\) 179 (-) Viral replication factor The T7 burst size is a good average as T7 does not have an unusually large or small burst size. [3]
\(\gamma\) 20 (1/h) Phage degradation rate - Based on the experimental results
\(\eta_c\) 0.0115 (1/h) Cry7Ca1 degradation rate Degradation of another cry toxin at pH 7.49 at 35 degrees [6]
\(\eta_R\) 2.1 (1/h) shRNA degradation rate - [10]
\(\alpha_c\) 1552 (-) Cry7Ca1 production per cell This value is based on the expression of the codon optimized sequence for Cry7Ca1, which replaced the gene 4.3 on the genome of the T7 phage. Pinetree simulation from our partnership with UT Austin
\(\alpha_R\) 1552 (-) shRNA production per cell Assumption that the same amount of shRNAs are produced as Cry7Ca1. Pinetree simulation from our partnership with UT Austin
\(X_{s,0}\) \(3.34 \cdot 10^{9}\) (number/mL) Initial concentration of susceptible cells Phocus is assumed to be able to infect all bacteria in gut -
\(X_{i,0}\) 0 (number/mL) Initial concentration of infected cells - -
\(P_0\) \(3.0 \cdot 10^{5}\) (number/mL) Initial phage concentration - -
\(T_0\) 0 (number/mL) Initial Toxin concentration - -
Details of shRNA calculations

According to literature, as little as 1 ng of dsRNA in the hemolymph is enough to produce a significant knockdown of locust genes [12]. Based on the length of the used dsRNA (447 bp) and the fact that dicer processes dsRNA into fragments of 20 bp, we estimated that this 1 ng of dsRNA corresponds to \(4.6 \cdot 10^{10}\) molecules of shRNA. As no literature on the volume of the desert locust midgut exists, we estimate it to be 0.13 mL (average length adult locust = 75 mm (A. Showler, 2008), width and height 5 mm (own estimate), ratio midgut:body = 0.08:1 (own estimate)). The model predicted a maximum shRNA concentration of \(1.21 \cdot 10^{12}\) molecules/mL, which then corresponds to a total of \(1.57 \cdot 10^{11}\) shRNAs molecules.

However, for shRNA to have an effect it first has to cross the gut lining and enter the hemolymph. In literature, we found that on average 1% of the initial shRNA concentration is able to cross the gut lining to enter the hemolymph within 5 minutes. After 30 minutes, this is around 4% of the initial value [10].

Thus, if \(1.57 \cdot 10^{11}\) molecules of shRNA are produced in the gut, approximately \(1.57 \cdot 10^9\) to \(6.29 \cdot 10^9\) molecules (1% - 4%) of shRNA should enter the hemolymph within 5 to 30 minutes of production.


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.

Details on LSA

A linear stability analysis is about characterizing the equilibrium points of a system. Therefore, determining the equilibrium points of the system is the first step of the linear stability analysis. An equilibrium point is defined as a point where all changes are zero and can therefore be found by setting the differential terms of the system (Eq. 1, Eq. 2, Eq. 3 and Eq. 4) to zero and solving the resulting system for the four variables, \(X_s\), \(X_i\), \(P\) and \(T\). Two equilibrium points (\(E_{-}\) and \(E_{+}\)) were found by hand as they are trivial solutions. The third equilibrium point (\(E_{*}\)) was more difficult to find by hand and we therefore used the SymPy (symbolic Python) solver (see our Python code).

\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}

The first two terms of the Taylor series of the model equations can be used as a linear approximation of the nonlinear system. The first term of the Taylor series just equals the function value of the point where the Taylor series is evaluated, which is zero for the equilibrium points of that system. Essentially, this only leaves the second term of the Taylor series (first-order derivatives). Therefore, the linear stability of each equilibrium point can be assessed by evaluating the equilibrium point in the Jacobian of the nonlinear system.

\begin{equation} J(X_s, X_i, P, T) = \begin{bmatrix}- P \epsilon + \mu \left(1 - \frac{X_{i} + X_{s}}{C}\right) - \frac{X_{s} \mu}{C} & - \frac{X_{s} \mu}{C} & - X_{s} \epsilon & 0 \\ P \epsilon & - \lambda & X_{s} \epsilon & 0 \\- P \epsilon & \beta \lambda & - X_{s} \epsilon - \gamma & 0 \\0 & \alpha \lambda & 0 & - \eta\end{bmatrix} \end{equation}

First, the Jacobian was evaluated in the \(E_{-}\) equilibrium point



\begin{equation} J(0,0,0,0) = \left[\begin{matrix}\mu & 0 & 0 & 0\\0 & - \lambda & 0 & 0\\0 & \beta \lambda & - \gamma & 0\\0 & \alpha \lambda & 0 & - \eta\end{matrix}\right] \end{equation}

The eigenvalues of this matrix are x1 = \(\mu\), x2 = \(-\lambda\), x3 = \(-\gamma\) and x4 = \(-\eta\). Since all involved parameters (growth rate, lysis rate, phage degradation rate and toxin degradation rate respectively) must be greater than zero to make sense biologically, the equilibrium point \(E_{-}\) is an linearly unstable equilibrium point.


Secondly, the Jacobian was evaluated in the \(E_{+}\) equilibrium point

\begin{equation} J(C,0,0,0) = \left[\begin{matrix}- \mu & - \mu & - C \varepsilon & 0\\0 & - \lambda & C \varepsilon & 0\\0 & \beta \lambda & - C \varepsilon - \gamma & 0\\0 & \alpha \lambda & 0 & - \eta\end{matrix}\right] \end{equation}

Solving the equation \(|J - xI| = 0\) to calculate the eigenvalues results in the following characteristic polynomial:

\begin{equation} (-\mu -x)(-\eta -x)(x^2 + (\gamma +C\varepsilon +\lambda)x + \lambda (C\varepsilon(1-\beta) + \gamma )) \end{equation}

This polynomial immediately shows that x1 = \(-\mu\) and x2 = \(-\eta\) are eigenvalues of the Jacobian evaluated in \(E_{+}\). The other two eigenvalues can be obtained by solving the remaining second order polynomial, which can be done using the ABC-formula:

\begin{equation} x_{3,4} = \frac{-(\gamma + C\varepsilon + \lambda) \pm \sqrt{(\gamma + C\varepsilon + \lambda)^2 - 4\lambda(C\varepsilon(1-\beta)+\gamma)}}{2} \end{equation}

As expected there are two solutions to this second order polynomial. However, when

\begin{equation} 4\lambda(C\varepsilon(1+\beta)+\gamma) = 0 \end{equation}

One of the eigenvalues vanishes (becomes zero). This expression becomes zero when

\begin{equation} \lambda = 0 \end{equation}

or when

\begin{equation} \beta = 1 + \frac{\gamma}{C\varepsilon} = \beta^* \end{equation}

In other words, only when there is no lysis or when \(\beta\) equals \(\beta^*\) then there are only 3 eigenvalues, which are all negative since all parameters have positive values:

\begin{equation} x_{3} = -(\gamma + C\varepsilon + \lambda) \end{equation}

which means that the equilibrium point \(E_{+}\) is linearly stable under these conditions. Furthermore, when:

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

then

\begin{equation} |-(\gamma + C\varepsilon + \lambda)| > \sqrt{(\gamma + C\varepsilon + \lambda)^2 - 4\lambda(C\varepsilon(1-b)+\gamma)} \end{equation}

which means that both eigenvalues \(x_3\) and \(x_4\) become negative. Thus, when \(\beta\) meets this condition then all eigenvalues are negative and hence \(E_{+}\) in linearly stable. However, when

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

then

\begin{equation} |-(\gamma + C\varepsilon + \lambda)| \lt \sqrt{(\gamma + C\varepsilon + \lambda)^2 - 4\lambda(C\varepsilon(1-\beta)+\gamma)} \end{equation}

which results in one eigenvalue being positive and one negative. Therefore, when \(\beta\) meets this condition, at least one eigenvalue is positive and hence \(E_{+}\) is linearly unstable.


At last, the Jacobian was evaluated in the equilibrium point \(E_{*}\):

\begin{equation} J(E_*) = \left[\begin{matrix}- \frac{\gamma \mu}{C \epsilon \left(\beta - 1\right)} & - \frac{\gamma \mu}{C \epsilon \left(\beta - 1\right)} & - \frac{\gamma}{\beta - 1} & 0\\\frac{\lambda \mu \left(C \beta \epsilon - C \epsilon - \gamma\right)}{C \beta \epsilon \lambda - C \epsilon \lambda + \gamma \mu} & - \lambda & \frac{\gamma}{\beta - 1} & 0\\\frac{\lambda \mu \left(- C \beta \epsilon + C \epsilon + \gamma\right)}{C \beta \epsilon \lambda - C \epsilon \lambda + \gamma \mu} & \beta \lambda & - \frac{\beta \gamma}{\beta - 1} & 0\\0 & \alpha \lambda & 0 & - \eta\end{matrix}\right] \end{equation}

The characteristic polynomial of this matrix will lead to a fourth order polynomial with very big coefficients, which makes calculating the Eigenvalues difficult. However, to be able to conclude on the linear stability of \(E_{*}\) different mathematical tricks exist, such as the Routh-Hurwitz criterion. The Routh–Hurwitz stability criterion is a mathematical test to determine whether all roots of the characteristic polynomial of a linear system have negative real parts. This analysis shows that when

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

the equilibrium point \(E_{*}\) is feasible and \(E_{+}\) becomes unstable. Furthermore, for increasing values of \(\beta\), the \(E_{*}\) equilibrium bifurcates towards a periodic solution. To mathematical proofs supporting these statements are rather extensive and therefore refer interested readers to the original paper of Beretta and Kuang [1].


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.

Detailed analysis PSA

Sensitivity of Cry7Ca1

We performed a PSA by ranging all parameters listed in Table 2, but also the initial susceptible bacteria concentration, from 20% to 500% one parameter at a time, determining its impact on the maximum Cry7Ca1 concentration (\(T_{C,max}\)) (Figure 5) (see our Python code).

Figure 4. 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 used parameter values for this analysis are found in Table 2.

Figure 4b shows that \(T_{C,max}\) is significantly impacted by changes in the amount of Cry7Ca1 toxin produced per cell (\(\alpha_C\)) and the initial concentration of susceptible cells (\(X_s\)). Moreover, the \(T_{C,max}\) seems to be almost linear dependent on both \(\alpha_C\) and \(X_s\). The \(T_{C,max}\) seems also to be affected by the Cry7Ca1 toxin degradation rate (\(\eta_c\)), lysis rate (\(\lambda\)) and carrying capacity (C) (Figure 4a), but not in the same order of magnitude as \(\alpha_C\) and \(X_s\). For \(\eta_c\), it is to be noted that the rate we used was based on the degradation rate of another Cry toxin, which was found in literature and measured in the soil at pH 7.49 and 35 ͦ C [6]. In contrast to the desert locust gut, the soil contains a limited amount of proteases [14]. Thus, in reality, the possibility exists that \(\eta_c\) exceeds this range, which may affect this prediction and the actual value for \(\eta_c\) could even be outside the 20-500% range of this sensitivity analysis. However, Cry toxins have evolved to work in the gut and, therefore, there are some natural mechanisms that protect them [7]. With this, we suggest that improving the stability of the Cry7Ca1 toxin, with for example protein engineering, would be an effective strategy for increasing the \(T_{C,max}\). However, future experiments should determine the actual stability of Cry7Ca1 in the locust gut to ensure the Cry7Ca1 degradation will not pose a problem to the effectiveness of PHOCUS. Increasing the (\(\lambda\)) also increases the \(T_{C,max}\) but leads quickly to an asymptotic plateau where increasing (\(\lambda\)) no longer increases the \(T_{C,max}\). This can be explained by the fact that the rate of infection (\(\varepsilon\)) limits the number of infected cells at a certain time and lysis can therefore not happen. Increasing the lysis rate would therefore not be the most effective strategy for increasing the \(T_{C,max}\). All other parameters seem to have a negligible effect on the \(T_{C,max}\).


Subsequently, as PHOCUS has to be fast, the sensitivity of the time to reach the maximum toxin concentration to all parameters was evaluated (Figure 5).

Figure 5. The time to reach the maximum concentration of Cry7Ca1 in the model of host-phage infection and toxin production and as a function of changes in parameters ranging from 20% to 500% of the original values. They were determined by running the system of ordinary differential equations for 100 hours. The used parameter values for this analysis are found in Table 2.

Figure 5b shows that the time to reach \(T_{C,max}\) is affected by \(\eta_c\) and \(\lambda\). As they increase, the \(T_{C,max}\) is reached faster. A larger \(\lambda\) causes the phages to lyse and infect more cells more rapidly, producing toxins faster as a result. Thus, increasing \(\eta_c\) could be a valid strategy to improve the speed of toxin production. However, in reality, changes in lysis rate are inversely correlated to changes in the viral replication factor and most likely also the number of toxins produced per cell [9]. Therefore, increasing the lysis rate could potentially also decrease the \(T_{C,max}\)and change the effectiveness of PHOCUS. Additionally, our model suggests that the toxin production already happens relatively fast. Thus, a strategy to increase the speed of PHOCUS would not be necessary.

Furthermore, Figure 5b shows that a larger \(\eta_c\) also decreases the time it takes to reach \(T_{C,max}\). However, higher degradation also slightly lowers the \(T_{C,max}\). Therefore, increasing the degradation rate does not necessarily lead to a faster toxin production.


Sensitivity of shRNA

The sensitivity of the maximum shRNA concentration (\(T_{R,max}\)) to parameter changes was also determined by individually ranging all parameters from 20% to 500% of the initial parameter value (Table 2) (Figure 6)

Figure 6. The maximum shRNA toxin concentration in the model of host-phage infection and toxin production and bacterial resistance as a function of changes in parameters ranging from 20% to 500% of the original values. The maximum toxin concentration was determined by running the system of ordinary differential equations for 100 hours. The used parameter values for this analysis are found in Table 2.

Similar to the PSA of Cry7Ca1, these results indicate that \(T_{R,max}\) linearly depends on the initial concentration of susceptible bacteria (\(X_s\)) and the amount of shRNA produced per bacterium (\(\alpha_R\)). In contrast to the PSA of Cry7Ca1, for the shRNA parameters the shRNA degradation rate (\(\eta_R\)) and lysis rate (\(\lambda\)) also seem to have a significant impact on \(T_{R,max}\). Because RNA is degraded much faster than the Cry7Ca1 toxin, increasing the stability of our shRNA could greatly increase \(T_{R,max}\) in the gut (Table 2). In our design, we therefore included different strategies to increase the stability of our shRNA (see Design). The effect of \(\lambda\) on \(T_{R,max}\) can be explained by the fact that with an increasing lysis rate, more shRNA can accumulate in the system. However, as mentioned before, increasing \(\lambda\) could lead to a decrease in \(\alpha_R\) as there is less time to actually produce the toxins. Moreover, increasing \(\alpha_R\) seems to have a bigger impact on \(T_{R,max}\) than increasing \(\lambda\). Therefore, increasing \(\alpha_R\) would be the preferred strategy to increase the effectiveness of our biopesticide rather dan increasing \(\lambda\).


In terms of time needed to reach \(T_{R,max}\), again \(\lambda\) and \(\eta_R\) have the most dominant impact, similar to the effect on Cry7Ca1 (Figure 7). However, our model shows that the speed of production is already relatively fast. Therefore, we argue that there is no need to further increase this speed and that focusing on increasing \(T_{R,max}\) is most important.

Figure 7. The time to reach the maximum concentration of shRNA in the model of host-phage infection and toxin production and bacterial resistance as a function of changes in parameters ranging from 20% to 500% of the original values. They were determined by running the system of ordinary differential equations for 100 hours. The used parameter values for this analysis are found in Table 2.

Here we showed that we need to maximise the number of cells susceptible to PHOCUS (by using a cocktail of pages) and the toxin number produced per cell (choosing a better insert site) to increase the maximum Cry7Ca1 concentration (\(T_{C,max}\)) . Furthermore, increasing the stability of Cry7Ca1 appears to have less of an impact on increasing \(T_{C,max}\). The lysis rate seemed to have the greatest impact on the time it takes to produce \(T_{C,max}\). However, our model suggests that the toxin production already happens relatively fast and therefore, a strategy to increase the speed of toxin production would not be necessary.

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.

Detailed explanation of model with resistance
Resistant bacteria population

All the differential equations and assumptions defined in section 1 of the model apply to this model. Additional assumptions were made for the description of phage resistance. First of all, the rate of change of resistance bacteria is assumed to increase as a result of logistic growth (first term), limited by the maximum population size (\(X_t = X_s + X_i + X_r\)) and evolution (second term). As mentioned before, the number of mechanisms by which bacteria can become resistant to a phage is endless. For simplification, we chose to coin these mechanisms in parameter \(\phi\), the rate of a bacteria acquiring resistance against a phage.

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

The addition of resistant bacteria to the model could change the equilibrium points of the system as well as their nature. Therefore, the new system was analysed again using a linear stability analysis. In total, the system contains two equilibrium points:

\begin{equation} E_{-} = (0,0,0,0, 0) \end{equation}
\begin{equation} E_{+} = (0, 0, 0, 0, C) \end{equation}

The Jacobian of this system is given by

\begin{equation} J = \left[\begin{matrix}- P \varepsilon + \mu \left(1 - \frac{X_{i} + X_{s} + X_{r}}{C}\right) - \frac{X_{s} \mu}{C} & - \frac{X_{s} \mu}{C} & - X_{s} \varepsilon & 0 & 0\\P \varepsilon & - \lambda & X_{s} \varepsilon & 0 & 0\\- P \varepsilon & \beta \lambda & - X_{s} \varepsilon - \gamma & 0 & 0\\0 & \alpha \lambda & 0 & - \eta & 0\\\phi - \frac{X_{r} \mu}{C} & - \frac{X_{r} \mu}{C} & 0 & 0 & \mu \left(1 - \frac{X_{i} + X_{r} + X_{s}}{C}\right) - \frac{X_{r} \mu}{C}\end{matrix}\right] \end{equation}

Subsequently, the Jacobian was evaluated in \(E_{-}\)

\begin{equation} J(0,0,0,0,0) = \left[\begin{matrix}\mu & 0 & 0 & 0 & 0\\0 & - \lambda & 0 & 0 & 0\\0 & \beta \lambda & - \gamma & 0 & 0\\0 & \alpha \lambda & 0 & - \eta & 0\\\phi & 0 & 0 & 0 & \mu \end{matrix}\right] \end{equation}

The eigenvalues of this matrix are x1 = \(\mu\), x2 = \(-\lambda\), x3 = \(-\gamma\), x4 = \(-\eta\) and x5 = \(\mu\). As some eigenvalues are positive, \(E_{-}\) is linearly unstable, just as in the system without resistance.

Subsequently, the Jacobian was evaluated in \(E_{-}\)

\begin{equation} J(0,0,0,0,C) = \left[\begin{matrix}0 & 0 & 0 & 0 & 0\\0 & - \lambda & 0 & 0 & 0\\0 & \beta \lambda & - \gamma & 0 & 0\\0 & \alpha \lambda & 0 & - \eta & 0\\- \mu + \phi & - \mu & 0 & 0 & - \mu\end{matrix}\right] \end{equation}

The eigenvalues of this matrix are all negative, which means that \(E_+\) is a linearly stable equilibrium point. These results show that, eventually, the system will converge to state where only resistance cells are present. For our biopesticide, the time scale at which the resistant cells take over the system is very important, because if this is too fast, our biopesticide will be less effective. Whether this is the case, depends on the numerical values for the parameters which are investigated in the next section.

Details on simulation with resistance

We solved the ODEs of the model numerically using the SciPy ordinary differential equation solver for Python (see our Python code). For most parameters, the same values as in Q1 (link) were used. Additional parameters and initial conditions can be found below (Table 3).

Table 3. Parameter and initial condition values for the prediction of the effect of bacterial resistance.
Parameter Value Explanation Source
\(\phi\) \(4.19 \cdot 10^{-9}\) (1/h) Bacterial resistance rate [17]
\(X_{s,0}\) \(3.34 \cdot 10^{9}\) (number/mL) Initial concentration of susceptible cells -
\(X_{i,0}\) 0 (number/mL) Initial concentration of infected cells -
\(P_0\) \(3.0 \cdot 10^{5}\) (number/mL) Initial phage concentration -
\(T_0\) 0 (number/mL) Initial Toxin concentration -
\(X_{r,0}\) 0 (number/mL) Initial concentration of resistant cells -

Continue to our 2D biofilm simulation!

References
  1. Beretta, E., & Kuang, Y. (1998). Modeling and analysis of a marine bacteriophage infection. Mathematical Biosciences. https://doi.org/10.1016/S0025-5564(97)10015-3
  2. Sezonov, G., Joseleau-Petit, D., &ampl; D’Ari, R. (2007). Escherichia coli physiology in Luria-Bertani broth. Journal of Bacteriology. https://doi.org/10.1128/JB.01368-07
  3. Nguyen, H. M., & Kang, C. (2014). Lysis delay and burst shrinkage of coliphage T7 by deletion of terminator Tφ reversed by deletion of early genes. Journal of virology, 88(4), 2107–2115. https://doi.org/10.1128/JVI.03274-13
  4. Abedon, S. T. (2016). Phage therapy dosing: The problem(s) with multiplicity of infection (MOI). Bacteriophage. https://doi.org/10.1080/21597081.2016.1220348
  5. Hunt, J., & Charnley, A. K. (1981). Abundance and distribution of the gut flora of the desert locust, Schistocerca gregaria. Journal of Invertebrate Pathology. https://doi.org/10.1016/0022-2011(81)90105-1
  6. Yao-yu, B., Ming-xing, J., & Jia-an, C. (2007). Impacts of Environmental Factors on Degradation of Cry 1 Ab Insecticidal Protein in Leaf-Blade Powders of Transgenic Bt Rice. Agricultural Sciences in China. https://doi.org/10.1016/S1671-2927(07)60031-5
  7. García-Gómez, B. I., Cano, S. N., Zagal, E. E., Dantán-Gonzalez, E., Bravo, A., & Soberón, M. (2019). Insect hsp90 chaperone assists bacillus thuringiensis cry toxicity by enhancing protoxin binding to the receptor and by protecting protoxin from gut protease degradation. MBio. https://doi.org/10.1128/mBio.02775-19
  8. Gibson B, Wilson DJ, Feil E, Eyre-Walker A. The distribution of bacterial doubling times in the wild. Proc Biol Sci. 2018 Jun 13 285(1880). pii: 20180789. doi: 10.1098/rspb.2018.0789
  9. Wang, I. N. (2006). Lysis timing and bacteriophage fitness. Genetics. https://doi.org/10.1534/genetics.105.045922
  10. Zhu, W., Vandingenen, A., Huybrechts, R., Vercammen, T., Baggerman, G., De Loof, A., Breuer, M. (2001). Proteolytic breakdown of the Neb-trypsin modulating oostatic factor (Neb-TMOF) in the hemolymph of different insects and its gut epithelial transport. Journal of Insect Physiology. https://doi.org/10.1016/S0022-1910(01)00086-5
  11. Wu, Y., Lei, C. F., Yi, D., Liu, P. M., & Gao, M. Y. (2011). Novel Bacillus thuringiensis δ-endotoxin active against Locusta migratoria manilensis. Applied and Environmental Microbiology. https://doi.org/10.1128/AEM.02462-10
  12. Wynant, N., Verlinden, H., Breugelmans, B., Simonet, G., & Vanden Broeck, J. (2012). Tissue-dependence and sensitivity of the systemic RNA interference response in the desert locust, Schistocerca gregaria. Insect Biochemistry and Molecular Biology. https://doi.org/10.1016/j.ibmb.2012.09.004
  13. Grosjean, H., & Fiers, W. (1982). Preferential codon usage in prokaryotic genes: the optimal codon-anticodon interaction energy and the selective codon usage in efficiently expressed genes. Gene. https://doi.org/10.1016/0378-1119(82)90157-3
  14. Vranova, V., Rejsek, K., & Formanek, P. (2013). Proteolytic activity in soil: A review. Applied Soil Ecology. https://doi.org/10.1016/j.apsoil.2013.04.003
  15. Song, H., Zhang, J., Li, D., Cooper, A. M. W., Silver, K., Li, T., & Zhang, J. (2017). A double-stranded RNA degrading enzyme reduces the efficiency of oral RNA interference in migratory locust. Insect Biochemistry and Molecular Biology, 86, 68–80. https://doi.org/10.1016/j.ibmb.2017.05.008
  16. Labrie, S. J., Samson, J. E., & Moineau, S. (2010). Bacteriophage resistance mechanisms. Nature Reviews Microbiology. https://doi.org/10.1038/nrmicro2315
  17. Capparelli, R., Parlato, M., Borriello, G., Salvatore, P., & Iannelli, D. (2007). Experimental phage therapy against Staphylococcus aureus in mice. Antimicrobial Agents and Chemotherapy. https://doi.org/10.1128/AAC.01513-06