2D biofilm simulation
The aim of PHOCUS is to develop a fast acting biopesticide. Using our first model, we showed that PHOCUS can produce high toxin concentrations within seven days. However, in that model we made the simplifying assumption that our system is always wellmixed, making it in effect ‘0dimensional’. In the desert locust gut, bacteria most likely reside in microcolonies [Dr. O. Lavy, personal interview, 1, 2]. The spatial organization of cells in these microcolonies can influence phage propagation and toxin production in the locust gut. We hypothesized that spatial effects, such as phage mobility and substrate diffusivity, may have an important influence on the performance of our biopesticide. We therefore modelled the temporal evolution of the number of bacteria, phages and produced toxin with a more detailed 2D biofilm simulation framework, to investigate the influence of these spatial effects on toxin production (Q1). Furthermore, as our biopesticide does not target all bacteria of the locust gut, we investigated the influence of biofilm heterogeneity on the resulting toxin production (Q2).
Model Explanation
We were inspired by a model simulation created by Simmons et al. [3]. We adapted this model to help us to answer our research questions (Figure 1).
In this model we simulate discrete bacteria and phages on a defined 2D grid. The growth of bacteria is governed by growthlimiting substrate that diffuses from the top of the 2D space into the biofilm, where substrate is depleted due to consumption by bacteria. Based on the substrate availability, bacteria grow and divide  filling up the grid. Once the number of bacteria on a given grid point reaches its maximum volume, the bacteria are redistributed to nearby grid points. Convection is modelled implicitly: the biofilm surface erodes in a heightdependent manner, reflecting the increase in shear force with distance from the surface (biomass erosion) and any free floating bacteria are removed from the system (biomass detachment). Phages move in the system according to a random walk, where the probability of moving depends on the interaction with the biofilm. Phages either associate with bacteria in the biofilm and initiate infection or diffuse into the bulk liquid, where they are removed from the system due to convection. Infected cells lyse after a defined time, releasing new phages and our locust specific toxins as they do so.
This biofilm framework simulates all these aforementioned events in discrete time steps (\(dt\)) along the following order:
 Diffusion of the nutrient substrate
 Biomass growth and division
 Lysis of infected bacteria and phage burst
 Erosion of biomass
 Phage movement
 Detachment of biomass
 Biofilm relaxation
 Phage infection
The toxins \(T\) were added to the system by counting the lysed cells (\(X_I\)) multiplied by the number of toxins produced per cell (\(\alpha\)) and are removed from the system with degradation rate \(\eta\) (Eq. 1).
As our toxins are small compared to bacterial cells and a single time step \(dt\) is relatively large (60 minutes), we assumed that the toxin diffusion is instantaneous compared to the other processes in the simulation. We therefore implemented the toxin production in a spaceindependent manner. This implementation should be sufficient to study the influence of spatial effects such as phage mobility (Q1) and biomass heterogeneity on the toxin production (Q2).
A more detailed explanation of the model can be found below:
Detailed model explanation
As this 2D biofilm simulation framework can be quite daunting, we aimed to provide an indepth explanation (with visual aid) of each step in the simulation cycle. In this way, we hope to give more insight into how such complex biofilm models work and how they could be implemented.
Diffusion of the nutrient substrate
The first step of the simulation is calculating the concentration profile of the substrate, which is described by the reactiondiffusion equation. The substrate diffuses with diffusion coefficient \(D_S\) and is consumed by bacteria according to a Monodtype substrate uptake reaction, where \(q_s\) is the maximum substrate uptake rate and \(K_S\) the half saturation constant (Eq. 2).
The nutrient concentration at the top of the system is kept constant (\(S_{max}\)) throughout the simulation (Dirichlet boundary condition). The boundary condition on the side and bottom of the system are equivalent to a zero gradient condition (Neuman boundary condition). It is assumed that the diffusion of substrate is instantaneous relative to bacterial and phagerelated processes and is therefore solved at a quasisteady state ( \( \frac{\partial S}{\partial t} = 0 \) ). In other words, a substrate gradient is computed on the 2D gridspace depending on the amount of bacteria in each grid (Figure 2). The substrate concentration per grid is then used in the subsequent iteration step to determine the increase in biomass (\(X\)).
Biomass growth and division
After determining the substrate concentration in each grid, the increase in biomass per bacteria is calculated. This is based on the amount of consumed substrate (\(q_s \cdot \frac{S}{S + K_S}\)) and the bacterial growth yield \(Y_{X/S}\).
The mass for each bacteria increases based on the change in biomass \( \frac{dX}{dt} \) and the time step \(dt\) of the iteration (Eq. 4) (Figure 3).
Once a bacterium reaches a critical mass \(m_{div} \), bacterial division occurs. The bacterial mass is divided in half and a second bacterium is added to that grid node (Figure 3).
Lysis of infected bacteria and phage burst
The third step of the simulation cycle determines which of the phage infected bacteria in each grid node burst due to cell lysis. Infected bacteria burst after a predefined lysis time \(\lambda \) and are removed from the grid node, releasing new phages. The number of new phages released depends on the burst size \( \beta \) (Figure 4). The new phages are free to diffuse throughout the biofilm and infect new cells.
Erosion of biomass
In reality, convection would result in shear stress and erode biomass on the surface of the biofilm. This is modelled implicitly by removing bacteria on the surface of the biofilm (Figure 5) with probability
where \(Er\) is the rate of erosion.
Phage movement
The fifth step of the simulation determines the phage movement through the biofilm. Phages move in the system according to a random walk, which is calculated as follows: First, the number of steps (\(n\)) that a phage could take in the next time interval (\(dt\)) is calculated by
where \( D_p \) is the diffusivity of the phage, \(dt\) is the time length of the iteration cycle and \(dl \) is the grid length scale. The time duration \( dt_p \) of each of those potential steps is then determined by
For each potential step (\(n\)), it is determined whether the phage
is removed from the system by convection.
To implicitly incorporate the loss of phages due to convection, any phages off the biofilm are removed with probability
\begin{equation} p = 1  e^{dt_p d^2 \delta_p} \tag{8} \end{equation}where \( d^2 \) is the distance away from the biofilm and \( \delta_p \) is the rate of removal due to convection.

remains in the current grid node or is be able to diffuse far enough to get into another grid node.
In 2 dimensions, the probability of diffusing into another grid node place is approximately 0.42, which is derived from the analytical solution of the diffusion equation in 2 dimensions (see [3] for derivation). The target node to which the phage could move is selected randomly.

has interacted with bacteria and therefore its motion ceases.
During each step that a phage takes, there is also a probability of the phage to interact with a bacterium to mimic the ability of the phage to bind to bacterial receptors, which inhibits the phage from further diffusion. The probability that the phage doesn’t interact, i.e. the probability that the phage is able to move to the target node, is described by
\begin{equation} p = 1  e^{dt_p (I_t + I_s)} \tag{9} \end{equation}where \( I_t \) is the impedance (measure for the amount of interaction between bacterium and phage) at the target grid node and \( I_s \) is the impedance at the source node. The impedance values are determined by
\begin{equation} I_x = \sum_i m^i_x I^i_x \tag{10} \end{equation}where \( m_x \) is the mass of an individual bacterium i and \( I_x \) is the impedance of that bacterium i with the phage. To put it simply, the more bacteria (which have a certain interaction with the phage) there are in a certain source and target grid node, the higher the probability that the phage interacts with one of those bacteria. Once the phage has interacted, its motion ceases.
As mentioned earlier, each iteration represents a discrete time step \( dt \). Phage mobility is thus described as the number of steps \( n \) allowed within \( dt \), wherein the duration of each step equals \( dt_p \). However, when a phage interacts with a bacterium, its motion ceases and a certain time \( dt_r \) remains. The remaining time \( dt_r \) is subsequently used to calculate the probability of infection (see phage infection section).
Detachment of biomass
The sixth step of the simulation involves the removal of biomass that is detached from the biofilm, which is another implicit implementation of convection. The detachment of bacteria can occur due to the creation of cavities where phages collect and propagate. Over several iterations the number of phages in these cavities increases and the number of bacteria that link the upper part of the biofilm with the main biofilm body decreases (Figure 7). These cavities keep growing until the link between the upper population and main biofilm body disappears, upon which the detached biomass is removed from the system.
Biofilm relaxation
Once a certain grid node reaches a maximum allowed biomass density \( X_{max} \), the biomass is redistributed to nearby grid nodes, which is the seventh step of the simulation. For each overfull gridnode, neighbouring grid nodes are explored until a suitable grid node is found with sufficient free volume. Subsequently, a line is computed towards this new grid node and the bacteria are redistributed along this line to satisfy the \( X_{max} \) condition.
Phage infection
The eighth step involves the infection of phages that interacted with a bacterium (see phage movement) with probability
where \( dt_r \) is the time remaining in the iteration after a phage has interacted with a bacterium and \( \gamma \) is the rate of infection. To put it simply, the more time there is for the phage to interact with the bacterium, the higher the probability of infection. Once infection occurs, the phage is removed from the grid node and the bacterium is now labelled as infected.
For more information about this 2D biofilm framework we refer readers to the original paper [3].
Q1. Which spatial effects are key for the performance of our biopesticide?
The ultimate goal of our biopesticide PHOCUS is to kill the locust within 7 days. To investigate the effect of spatial effects on our biopesticide, we analysed the influence of phage mobility on the maximum produced number of toxins toxin. For increasing values of the impedance \( I \), the ability of the phage to diffuse through the biofilm decreases. We performed a parameter sweep varying the impedance \( I \) and burst size \( \beta \) (Figure 10). The main result of this parameter sweep is that the burst size (\( \beta \)) could be used to compensate for impaired phage mobility.
Therefore, selecting a high burst size phage is an effective strategy to compensate for impaired phage mobility and thus results in a higher toxin production.
To further investigate which strategies could be used to compensate for impaired phage mobility, we performed a second parameter sweep to study if increasing the initial phage concentration would increase the maximum number of produced toxins (Figure 11).
This sweep again demonstrates the importance of phage mobility to toxin production and indicates that phage propagation and toxin production do not depend much on the initial phage concentration. This result indicates that biofilm structures actually promote phage propagation, as all bacteria are packed relatively close together. For our biopesticide, this implies that only a small number of phages would have to be delivered to a microcolony in the gut to induce toxin production. This makes PHOCUS feasibly suited for application using the UltraLowVolume standard of the FAO.
Q2. What is the influence of biofilm heterogeneity on the resulting toxin production?
Our bacteriophagebased biopesticide will most likely not target all bacteria of the locust gut. Our second goal of this model was to investigate the effect that this biomass heterogeneity could have on the toxin production. Therefore, we performed additional simulations where a second bacterial species was introduced in the system. This species had exactly the same bacterial properties as the first species, but that was nonsusceptible to the phage. During these simulations, we varied the fraction at which susceptible/nonsusceptible cells were introduced in the system. For each fraction, the simulation was repeated 3 times. Figure 12 shows the result of such a simulation that started with 5% susceptible and 95% nonsusceptible cells.
From this result, we could qualitatively observe that in the simulation, the susceptible cells (red) were not getting reached while there were still phages in the system. In this simulation, this led to almost no toxin production (Figure 13, yellow line), while the other two replicate simulations led to an expected toxin production profile (Figure 13, green and blue lines).
From this we can qualitatively observe that low fractions of susceptible bacteria in the biofilm leads to a decrease in the probability of phages to actually propagate before being removed from the system by convection. We only observed this qualitative behaviour for susceptible bacteria fractions smaller than 10%. For higher fractions we did not observe any difficulty for the phage to propagate and produce toxins, which again demonstrates the strength of our biopesticide in toxin production in such biofilm structure. For our biopesticide this would imply that, if the fraction of bacteria we target is too low, then the effectiveness of our biopesticide can be reduced. PHOCUS should therefore target a minimum fraction of the bacteria in the gut to prevent loss of effectiveness.
Due to the probabilistic nature of some processes in this framework, we were only able to draw qualitative conclusions from our simulations. To strengthen the observed qualitative behaviour we suggest to increase the amount of simulation repeats to eliminate any probabilistic variance. Furthermore, to be able to more quantitatively conclude on the effect of the fraction of susceptible/nonsusceptible bacteria, we suggest performing an extensive parameter sweep varying this fraction. From this data, it can be determined if the relation between fraction of susceptible cells and the resulting produced toxin is linear or if spatial constraints prevent the phages from infecting the susceptible bacteria at lower fractions.
What did we learn?
With this model, we demonstrated the importance of phage mobility in the ability of our biopesticide to produce high toxin quantities. Interestingly, we showed that selecting a phage with a higher burst size could help overcome limitations in phage mobility and would therefore be a good strategy to increase toxin production. This observation differs from the results of our first model, where we showed that increasing the burst size has a negligible effect on the maximum produced number of toxins (see model 1). From this, we conclude that phage infection and toxin production in spatially structured biofilm contexts fundamentally differs from ideallymixed systems, as assumed in our first model. Furthermore, we showed that the initial amount of phages introduced into the system has a negligible effect on toxin production. From this, we conclude that biofilm structures promote phage propagation as cells are closely packed together. For PHOCUS, this is an advantage as only a small number of phages has to be delivered to a microcolony in the gut to induce toxin production.
Additionally, we showed that biomass heterogeneity, i.e. a biofilm consisting of both susceptible and nonsusceptible bacteria, can lead to a decrease in the probability of phages to actually propagate before being removed from the system by convection. Our biopesticide should therefore target a minimum fraction of the bacteria in the gut to prevent loss of effectiveness.
Simulation details
Simulation initialisation and parameters
Each simulation was initialised on a grid of 450 \(\mu m \) (xaxis) by 90 \(\mu m \) (xaxis), where each grid node had a size of 3 \(\mu m \) by 3 \(\mu m \). Within each grid, bacteria and phages are tracked individually but assumed to interact randomly. Each simulation was initialised by placing 150 bacteria randomly on the bottom of the grid, where they can grow and multiply based on the diffusion of substrate from the bulk layer above the biofilm. As our goal was to investigate the influence of spatial effects on phage infection and toxin production, we adapted the same parameters as used in the original biofilm simulation framework [3], which are listed below (Table 1):
Parameter  Value  Description  Source 

\(x_{max}\), \(y_{max}\)  450 \(\mu m \), 90 \(\mu m \)  Physical size of the system   
\({dl}\), \(dV\)  3 \(\mu m \), 27 \(\mu m^3 \)  Length and volume of a grid element   
\( S_{max} \)  \(6\) (mg/L)  Bulk substrate concentration  [3] 
\(D_S\)  \( 2.3 \cdot 10^{6} \) cm^{2}/s  Diffusivity of substrate  [3] 
h  15 \(\mu m \)  Diffusion boundary layer height  [3] 
\(K_N\)  \(1.18\) (mg/L)  Halfsaturation constant for substrate  [3] 
\(m_s\)  \( 10^{12} \) g  Bacterial mass per cell  [3] 
\(q_s\)  28.5 g/day  Substrate uptake rate  [3] 
\(X_{max}\)  200 g/L  Maximum biomass density  [3] 
\(Y_{X,S}\)  0.495  Biomass yield on substrate  [3] 
\(\beta\)  120 ()  Phage burst size  [3] 
\(D_P\)  \( 3.82 \cdot 10^{7} \) cm^{2}/s  Phage diffusivity  [3] 
\(I\)  \( 6 \) (g day)^{1}  Rate of interaction (impedance) of phage particles with biomass particles  [3] 
\(\delta_P\)  0.10 ( \(\mu m^2\) h)^{1}  Phage removal rate  [3] 
\(\lambda\)  28.8 min  Lysis time  [3] 
\(\gamma\)  2.92 h^{1}  Infection rate per biomass per phage  [3] 
\(dt\)  60 min  Timestep of simulation   
\(m_{div}\)  \( 1.33 \cdot 10^{12}\)  Bacterial division mass  [3] 
\(\eta\)  0.0115 (1/h)  Toxin degradation rate  [4] 
\(\alpha\)  1552 ()  Toxin production per cell  Pinetree simulation from our partnership with UT Austin 
Initialising phage infection
Before phages are added to the system, the biofilm is first given some time to grow. We initialised phage infection once the biofilm reached a height of 10 \(\mu m \). We based this estimate on a light micrograph image of a section through two rectal pads of the desert locust [2]. In each standard simulation, 450 phages were distributed randomly across the surface of the biofilm to initiate infection. Each simulation ended 4 days after phage initialisation.
Details on parameter sweeps
To also study the interplay between parameters, two parameters were varied throughout a sweep, which resulted in a heatmap as shown in the results section. In the first sweep, the phage burst size (30, 60, 120, 180, 240) and the impedance (3, 4.5, 6, 7.5, 9) were varied (Figure 10). In the second sweep, the initial phage number (10, 110, 220, 330, 640) and the impedance (3, 4.5, 6, 7.5, 9) were varied (Figure 12). The simulation for each parameter set was repeated 3 times (triplicate). Depending on the used parameters, the toxin profile could contain high fluctuations making it more difficult to obtain the maximum toxin value. We therefore used a mean sliding window to smoothen the toxin profiles and obtain the mean maximum toxin value for each unique parameter combination (Figure 14).
High performance computer cluster
The computational load of Model 2 is high as a single simulation requires 3060 minutes to execute. Performing a parameter sweep on such lengthy simulation would take a long time to execute. Fortunately, we were able to use one of the HighPerformance Computer Cluster at our university, allowing us to run many simulations simultaneously. We also wish to express our gratitude for the help provided by Emilia Simmons, as she provided us with simulation examples and a starter guide explaining how to execute this framework on a computer cluster using the distributed resource manager TORQUE. The (adjusted version) of the framework in Python, including the guide provided by Emilia various scripts for data analysis can be found here.
Continue to our phage infection and toxin production model!
References
 Dillon, R., & Charnley, K. (2002). Mutualism between the desert locust Schistocerca gregaria and its gut microbiota. In Research in Microbiology (Vol. 153, Issue 8, pp. 503–509). Elsevier Masson. https://doi.org/10.1016/S09232508(02)01361X
 Hunt, J., & Charnley, A. K. (1981). Abundance and distribution of the gut flora of the desert locust, Schistocerca gregaria. Journal of Invertebrate Pathology, 38(3), 378–385. https://doi.org/10.1016/00222011(81)901051
 Simmons, M., Drescher, K., Nadell, C. D., & Bucci, V. (2018). Phage mobility is a core determinant of phagebacteria coexistence in biofilms. ISME Journal, 12(2), 532–543. https://doi.org/10.1038/ismej.2017.190
 Yaoyu, B., Mingxing, J., & Jiaan, C. (2007). Impacts of Environmental Factors on Degradation of Cry 1 Ab Insecticidal Protein in LeafBlade Powders of Transgenic Bt Rice. Agricultural Sciences in China. https://doi.org/10.1016/S16712927(07)600315