Team:TUDelft/Model/Biofilm Simulation

PHOCUS

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 well-mixed, making it in effect ‘0-dimensional’. 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).

Figure 1. Example output of biofilm simulation framework with phage infection created by Simmons et al. [3]. Bacteria (red) and phages (black) are modelled discreetly on a 2D grid. After a defined period of bacterial growth on growth-limiting substrate, phages are introduced in the system which can infect and lyse bacteria (green), releasing new phages and toxins as they do so.

In this model we simulate discrete bacteria and phages on a defined 2D grid. The growth of bacteria is governed by growth-limiting 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 height-dependent 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:

  1. Diffusion of the nutrient substrate
  2. Biomass growth and division
  3. Lysis of infected bacteria and phage burst
  4. Erosion of biomass
  5. Phage movement
  6. Detachment of biomass
  7. Biofilm relaxation
  8. 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).

\begin{equation} \frac{dT}{dt} = \alpha X_{lysed}(t) - \eta T \tag{1} \end{equation}

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 space-independent 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 in-depth 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 reaction-diffusion equation. The substrate diffuses with diffusion coefficient \(D_S\) and is consumed by bacteria according to a Monod-type substrate uptake reaction, where \(q_s\) is the maximum substrate uptake rate and \(K_S\) the half saturation constant (Eq. 2).

\begin{equation} \frac{\partial S}{\partial t} = D_S \nabla^2 S - \frac{q_s S}{S + K_S}X \tag{2} \end{equation}

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 phage-related processes and is therefore solved at a quasi-steady state ( \( \frac{\partial S}{\partial t} = 0 \) ). In other words, a substrate gradient is computed on the 2D grid-space 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\)).

Figure 2. First step of the simulation is determining the substrate (green circles) concentration profile by solving a Monod-type reaction-diffusion equation at quasi-steady state (\( D \nabla^2 S + \frac{q_sS}{S + K_S}S = 0 \)). The computed substrate gradient is used in the subsequent step of the simulation to determine the increase in biomass in each grid (brown squares).
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}\).

\begin{equation} \frac{dX}{dt} = Y_{X,S} \cdot q_s \cdot \frac{S}{S + K_S} \tag{3} \end{equation}

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

\begin{equation} X_{new} = X_{old} + \frac{dX}{dt} \cdot dt \tag{4} \end{equation}

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

Figure 3. Second step of the simulation is biomass growth, where a bacterium (light brown circle) inside a grid node (dark brown squares) grows depending on the available substrate S, substrate uptake rate \(q_s\), growth yield \(Y_{X/S}\) and timestep \(dt\) of the cycle. A bacterium divides when it reaches the critical mass \(m_{div} \).
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.

Figure 4. Third step of the simulation is lysis of infected bacteria (brown circle) and phage burst. Each infected bacterium burst after a predefined lysis time \(\lambda\), releasing new phages (green) into the system.
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

\begin{equation} p = 1 - e^{-E_r \cdot dt} \tag{5} \end{equation}

where \(Er\) is the rate of erosion.

Figure 5. Convection is modelled implicitly by removing bacteria (dark brown circles) on the biofilm surface with probability \(1 - e^{-E_r \cdot dt}\), where \(E_r\) is the rate of erosion and \(d_t \) the duration of the iteration.
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

\begin{equation} n = \frac{D_p dt}{2dl^2} \tag{6} \end{equation}

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

\begin{equation} dt_p = \frac{dt}{n} \tag{7} \end{equation}

For each potential step (\(n\)), it is determined whether the phage

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

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

  3. 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).

Figure 6. Fifth step of the simulation is calculating phage (green) movement based on a random walk. During this movement phages can be removed from due to convection when they are off the biofilm, diffuse to another grid node (brown squares) or interact with a susceptible bacterium (dark brown circles) and cease motion. Within the timestep \(dt\) phages can take a total of \(n\) steps. The duration of each steps equals \( dt_p \). When a phage interacts with a bacterium, its motion ceases and the remaining time of the iteration \(dt_r\) is used to calculate the probability of infection (see phage infection).
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.

Figure 7. The removal of biomass that is detached from the biofilm (brown squares), to implicitly simulate convection.
Biofilm relaxation

Once a certain grid node reaches a maximum allowed biomass density \( X_{max} \), the biomass is redistributed to near-by 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.

Figure 8. The eight step of the simulation is the redistribution of bacteria (circles) to neer-by grid nodes (squares), once a certain grid node reaches a maximum allowed biomass density \( X_{max} \). For each overfull grid node (red), neighbouring grid nodes are explored until a suitable grid node is found with sufficient free volume (green). 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

\begin{equation} p = 1 - e^{-dt_r \gamma} \end{equation}

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.

Figure 9. The ninth step of the iteration involves determining whether a phage (green) that interacted with a bacterium (brown circle) is able to infect with probability \(p = 1 - e^{-dt_r \gamma} \). \( dt_r \) is the time remaining in the iteration with length \( d_t \) after a phage interacted with a bacteria and \( \gamma \) is the rate of infection.

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.

Figure 10. Heat map of a parameter sweep over impedance (x-axis) and burst size (y-axis). The values in each grid correspond with the mean maximum produced number of toxins toxins of 3 repeated simulations. When impedance decreases, the phages can more easily diffuse through the biofilm. This increases the speed of phage propagation and results in a higher toxin production. Similarly, increasing the phage burst size speeds up phage propagation and results in a higher toxin production.

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

Figure 11. Heat map of a parameter sweep over impedance (x-axis) and initial phage concentration (y-axis). The values in each grid correspond with the mean maximum number of produced toxins of 3 repeated simulations. Increasing the initial number of phages to induce infection has a negligible effect on the mean maximum number of produced toxins. This result demonstrates that the spatial organization of cells in biofilm structures promotes toxin production, as phages can easily spread when bacteria are packed closely.

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 Ultra-Low-Volume standard of the FAO.


Q2. What is the influence of biofilm heterogeneity on the resulting toxin production?

Our bacteriophage-based 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 non-susceptible to the phage. During these simulations, we varied the fraction at which susceptible/non-susceptible 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% non-susceptible cells.

Figure 12. Biofilm simulation initialised with 5% susceptible (red) and 95% non-susceptible cells (blue). X-axis and y-axis both represent the physical space of the simulation in μm. A) 450 phages (black) were introduced in the system after 1.71 days of bacterial growth. B) After less than a day, all phages where removed from the system due to convection (number of phages = 0), without the occurrence of phage infection. This shows that for low fractions of susceptible cells, the probability of phages to initialize infection decreases, which could negatively affect the performance of our biopesticide.

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

Figure 13. Toxin profiles of 3 biofilm simulations (triplicate) initialised with 5% susceptible and 95% non-susceptible cells. Two of the replicates showed an expected toxin production production profile in time (blue and green). During the third simulation, however, no toxins were produced (yellow). This shows that for low fractions of susceptible cells, the probability of phages to initialize infection decreases, which could result in no toxin production.

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/non-susceptible 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 ideally-mixed 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 non-susceptible 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 \) (x-axis) by 90 \(\mu m \) (x-axis), 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):

Table 1. Parameters and initial conditions for model prediction.
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} \) cm2/s Diffusivity of substrate [3]
h 15 \(\mu m \) Diffusion boundary layer height [3]
\(K_N\) \(1.18\) (mg/L) Half-saturation 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} \) cm2/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).

Figure 14. We used a mean sliding window using 15 data points to obtain the mean toxin value of 3 simulations of the same unique parameter combination.

High performance computer cluster

The computational load of Model 2 is high as a single simulation requires 30-60 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 High-Performance 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
  1. 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/S0923-2508(02)01361-X
  2. 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/0022-2011(81)90105-1
  3. Simmons, M., Drescher, K., Nadell, C. D., & Bucci, V. (2018). Phage mobility is a core determinant of phage-bacteria coexistence in biofilms. ISME Journal, 12(2), 532–543. https://doi.org/10.1038/ismej.2017.190
  4. 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