Team:William and Mary/Model
Mathematical Modeling
Abstract
We constructed mathematical models that explored three different aspects of biofilms: biofilm formation, diffusion of chemical signals, and Turning pattern formation within biofilms. By generating, testing, and fitting parameters using our differential equations, we were able to provide the wet lab team with information essential for their experiments including but not limited to, quantitating percentage of adherent cells, calculating the concentration of chemical signals in certain areas, estimating dynamics of chemical signals such as production, degradation and diffusion in our system.
What We Did
- We developed our first set of models based on the ODE equation systems and parameters from Jin and Riedel’s paper Biofilm Lithography Enables High-resolution Cell Patterning Via Optogenetic Adhesin Expression. We chose to build these 2 models because of their abilities to provide explanations of the mechanism behind light-regulated biofilm formation and insights of related natural biofilm formation processes. We then refined and extended these models by adding in factors such as cell collision, cell growth, and cell death.
- We constructed a diffusion model with ordinary differential equations to predict the interactions results between two adjacent biofilms via chemical substance production and diffusion. This model is an example of distance-dependent patterning. We scaled our simulations using the size of the biofilm that the wetlab team produced. We also modeled various concentration of IPTG and various length of IPTG incubation to help the wetlab team revise their future double ring experiments.
- We combined our diffusion model with Turing equations to explore the possibilities of forming Turing patterns within our biofilm diffusion system. We have readjusted our ordinary differential equations multiple times in order to better simulate the pattern the Turing patterning would produce in nature, but we did not quite capture the intricate relationship between the fast diffusing inhibitor and the slow diffusing activator within our model. However, we can observe from the results that with each modification of the ordinary differential equations, our results were becoming more stochastic and closer to what actual Turing pattern would look like.
Motivations
We chose to model light-regulated biofilm formation because it gave us a comprehensive understanding of the mechanism behind biofilm formation. By utilizing this mechanism, we could precisely control the growth of biofilms in both spatial and temporal dimensions. We also want to test their hypothesis regarding light-regulated biofilm formation before we start our experiment. They hypothesized that the direct cause MG1655/pDawn-Ag43 biofilm patterning was a decrease in their desorption/detachment rate, in contrast to an increase in their adsorption/attachment rate.
The models are designed to simulate how a given number of cells would behave, in terms of attachment and detachment to the surface, when the only variable in the environment is the presence of blue light. With the ODE equations and parameters they provided in their paper, we further improved the accuracy of these models by adding in factors such as cell growth and cell death. This has made this simulation much more comparable to real-life experiments. With the success of our simulations, we can grasp the theories behind the process of biofilm-forming under blue light, which gives us a firm foundation to proceed to build our chemical diffusion model.
Our models are potentially applicable for studies in natural microbial communities as well as studies in engineering living biomaterials.
Models / Simulations
We used the Monte-Carlo method to simulate cell swimming and adsorption/desorption. We initialized a group of individual cells, assigned them properties such as adhesin level, velocity, direction, position on the two-dimension plane, and whether or not they were illuminated by light based on their position. We allowed those individual cells to randomly move, turn, desorb/detach, and adsorb/attach in a virtual chamber with one side illuminated by blue light and the other side left dark. The movements and turnings of the cell were random and were the same for every single cell, but the attachment and detachment of cells were regulated by ordinary differential equations. After a given period, our models would display positions that have cells attached to it as a dot. By observing how many dots there were in a region, we were able to tell if there is a biofilm formed in this region.
To test their hypothesis, we created 2 different models. Model 1 corresponds to the rejected hypothesis, in which sets the detachment rate to be constant and attachment rate to be an increasing function of Adhesin. Model 2 corresponds to the verified hypothesis, in which sets the attachment rate to be a constant, and the detachment rate to be a decreasing function of Adhesin. Later, our team deemed this model to be incomplete because it did not take account of cell collision, cell division and cell death. We fixed these issues by reorienting the cell's directions and readjust cells’ positions on the surface when two or more cells have overlapping positions. We modeled cell division and cell death by periodically remove a random group of cells and creating a group of new cells entries, the new cells all have an initial adhesin level of 1 because that is the starting Adhesin level we have assigned to all cells at the beginning of the simulation.
Model 1 and 2 are implemented in Matlab using a Forward Euler’s numerical approach, with the time being discretized to dt = 0.1s timesteps. Simulations were run for 16 hours over an area of 600µm². Exactly like the numbers that Jin and Riedel chose to use. In the first step of each model, the adhesin level for each individual cell is regulated and updated in each timestep by an ordinary differential equation:
The second step of each model consists of updating the positions of cells based on the current position, velocity and direction based on a combination of tumbling probability and tumbling frequency. If the simulation has determined, based on the combination, that the cell should move in this timestep, the cell will be given a random direction, and the position of the cell will be updated using the equations below:
We set the velocity of cells to be constant in here with a numerical value of 14 µm/s. The direction of cells are randomly generated numbers within the range of (0, 2*pi)
The third and final step of model 1 and model 2 are different. In our simulations, we used blue light to illuminate the area with the abscissa from 0 to 400 and the ordinate from 0 to 600. In model 1, where the cell attachment rate is increased upon illumination, the relationship were set as follows:
In model 2, where the cell detachment rate is decreased upon illumination, the relationship were set as follows:
The detailed explanation of each variable and their corresponding numerical value could be found in the table below
After this iteration was done, one timestep of the simulation would be finished. In our simulation, we chose to simulate 16 hours in terms of experiment times, which means the above three steps would be repeated 576,000 times.
Results
The final result that we got for model 1 is shown below:
The final result we got for model 2 is shown below:
The result of those two models confirmed the hypothesis that Jin and Riedel have made regarding light-regulated biofilm formation: The pDawn-Ag43 mediated patterning occurs because of its significant impact on decreasing the detachment rate of the already attached cells.
Conclusion
The successful results from this model have demonstrated to us the mechanism behind light-regulated biofilm formation. The results are also very useful for us in the future because we can engineer biofilms with the shapes and sizes that we desire. The result of this modeling is for the diffusion model and our explorational model of Turing equations.
How it influenced wetlab: by successfully recreating the biofilm lithography model, we helped the wet lab to verify the validity of this method of biofilm patterning. The biofilm later confirmed our model by conducting experiments and showing results that corresponds to our modeling results.
Reference
Jin, X., & Riedel-Kruse, I. H. (2018, April 3). Biofilm Lithography enables high-resolution cell patterning via optogenetic adhesin expression.
Motivations
We drew our inspiration for this model from researcher Pulin Li’s paper “Morphogen gradient reconstitution reveals Hedgehog pathway design principles”. We chose to model chemical substance diffusion within biofilm because we want to use synthetic biology techniques to investigate the cell-cell communication systems, thereby allowing us to explore novel forms of patterning. This form of patterning is important to us because it is a form of self-patterning by the biofilm itself. If we could engineer the biofilm to constantly induce its outer areas to produce more biofilm, then the result would be a self-organize, stable biofilm that could be demonstrated useful for many areas.
Our model is designed to simulate the diffusion and transmission of a certain chemical signal between two different strains of E. coli. In our model, we first utilized the biofilm lithography method described in section 1 to create a double ring form of biofilm, consisting of two different strains of E. coli. A graphical example is shown below to illustrate the biofilm lithography setup better.
As you can see in this double ring image. Both biofilms were formed using the biofilm lithography method. With the center green dots representing the sender strain, and the outer yellow dots representing the receiver strain. It is worth noting that in our experiment, the sender strain was placed down first, and then all the peripheral sender strains were gently washed off. We tried to capture this process in our model by removing cells that were not in the designated biofilm formation region. This explains the crisp clear boundary between the sender strain and the receiver strain.
Model / Simulations
The biology behind our model is as follows: In this model, the sender strain, denoted in the graph by the color green, when chemically induced by IPTG, will produce the signal homoserine lactone and diffuse outwards. The receiver strain, denoted in the graph by the color yellow, will respond to this HSL by inducing the production of mScarlet. Thereby allowing us to observe the diffusion result of HSL directly.
The chemical diffusion model employed the Forward Euler’s numerical approach, with the time being discretized to dt = 0.1s timesteps. This simulation was run for 16 hours over an area of 6000(micrometer)^2. We chose those numbers because it is representative of the real biofilm size from our experiment. We first initialized the individual cells for strain 1 biofilm. For every single cell that belongs to strain 1, it still had its adhesin level, velocity, direction, and position. We decided to add an HSL concentration feature to each strain 1 cell to keep track of how much HSL that particular cell had produced over time. For every single strain 2 cells, it has its adhesin level, velocity, direction, position and its mScarlet level that keeps track of how much mScarlet it has produced over time.
For each time step, the model would go through two processes that update different values. During the first step, The concentration of HSL and mScarlet were updated using the following equations:
The explanations for each variable are provided in the table below:
The second step of the model consists of solving a diffusion equation. The steps below will walk you through our process of resolving a complicated equation utilizing a finite difference approach.
Often, the two-dimensional diffusion is characterized by Fick’s second law:
The explanation for each variable is provided in the table below:
Since we are using a forward Euler’s approach, the time has already been discretized, and we now need to discretize space so that this ordinary differential equation would become a system of linear equations. To do that, we divided the plane that our biofilms are located on into individual grids, as shown in the graph below.
Thus, we have discretized space. Our diffusion equation now could be simplified as
The table below provides an explanation for each variable:
Essentially, we summed up the concentration difference between our target grid and the 8 grids surrounding our target grid. This summation determined whether HSL should be diffused out of the target grid or diffused into the target grid. We decided to place the degradation term in this equation because degrading by position produces more accurate results than degrading by individual cells.
In summary, for each time step, we iterated through all the grids in this plane. For each grid, we would count the amount of HSL in it, calculate the sum of differences between itself and its surrounding grids, and make the decision of either diffusing HSL outwards or inwards using the simplified diffusion equation. We admit that it looks like a long run time, but it is a simple and powerful computational method that provides us accurate results.
Using this method, we collected the fluorescence intensity for the whole system after 5, 10 and 16 hours. The comparison of the HSL concentration level in those three periods are shown below.
5 hours
10 hours
16 hours
The comparison of the fluorescence intensity level in those three time periods are shown below.
5 hours
10 hours
16 hours
Conclusion
The result of this model showed us the successful diffusion of chemical substances and the induced fluorescence level as a response to that. This confirmed the hypothesis that our Wetlab has made that the fluorescence should be strongest in a ring surrounding the circular sender biofilm, and this ring of fluorescence should fade outwards in a gradient.
How it influenced wetlab: We have been communicating with wetlab about our modeling result back and forth. Our Model was initially in the range of micrometers using the parameters from Stanford researcher's group. However, after communicating our results with the wetlab team, we realized that the real biofilm size was at least 10 times larger than the scale that we are using right now. Thus we expanded our dimensions to adapt to their experimental conditions. We later helped the wet lab team with their experiment by modeling different IPTG input as a function of time and also various lengths of IPTG incubations. By observing the resulting fluorescence level, we gave the wetlab team a desirable range for choosing their parameters for future ring experiments.
Reference
Li, P., Markson, J. S., Wang, S., Chen, S., Vachharajani, V., & Elowitz, M. B. (2018, May 4). Morphogen gradient reconstitution reveals Hedgehog pathway design principles.
Motivation
For this model, we took inspiration from David Karig and Michael Martini’s paper Stochastic Turing patterns in a synthetic bacterial population.
We chose to make this exploration step of combining our diffusion model with Turing equations because we want to explore the possibility of design and control biofilms using the Stochastic Turing patterning method. Compared with other patterning methods such as printing and pattern substrates, Stochastic Turing patterns has the advantage of not being limited by external management such as access to substrate, expenses, or surface pretreatment that are often required for other patterning methods. The minimal demand for external interference in pattern formation has made Stochastic Turing pattern a very desirable target for us to model.
To incorporate the Turing patterns within our diffusion model, we took considerations of the circuit element that we already have and decided to modify the sender strain. The biofilm setup is the same as in the diffusion model, but the sender strain now is activated to produce both the activator and inhibitor signals simultaneously with IPTG input as a function of time. The following process happened as we ran the model: with the input of IPTG, the sender strain activated plasmid PFNK-512, which activated the activator and inhibitor signals. With the production of activator signals, PFNK-816, which inhibited the production of activator and inhibitor signals were activated and inhibited the production of both signals in the circuit. The activator and inhibitor signals were changed during each time-step depending on both the amount of PFNK-512 and the amount of PFNK-816 here. In our model, we simplified this circuit by making activator responsible for both its own and the inhibitor’s production, and we also made the inhibitor responsible for inhibiting both itself and the activator. The final result of this model is represented by a concentration mapping of fluorescent protein, whose production is only activated in the presence of the activator signals.
Our model is designed to simulate the diffusion, transmission and interaction of the activator and inhibitor chemical signals. In our model, we first utilized the biofilm lithography method described in section 1 to create a double ring form of biofilm, consisting of two different strains of E. coli. A graphical example is shown below to illustrate the biofilm lithography setup better.
Note that to emphasize the interaction and diffusion of those chemical substances. We chose to expand the strain 2 biofilms so that there were more space for chemical signals to diffuse and interact.
This model employs the Forward Euler’s numerical approach, with the time being discretized to dt = 0.1s timesteps. This simulation was run for 16 hours over an area of 6000(micrometer)^2. For every strain 1 bacterium, it has its adhesin level, velocity, direction, position, activator signal concentration, and inhibitor signal concentration. Whereas every single strain 2, bacterium has its adhesin level, velocity, direction, position, and fluorescent protein concentration.
For each time step, the model will go through two processes, the first step will update the concentration of activators and inhibitors for strain 1 bacteria, and the concentration of fluorescent protein for strain 2 bacteria. The concentration for activator, inhibitors, and DsRed(fluorescent protein) were updated using the following equations:
The explanations for each variable are provided in the table below:
In the second step, the model used the same finite difference method to resolve the diffusion process of activators and inhibitors. We employed the same grid method to keep track of the activator concentration and inhibitor concentration in every single location separately. Since we now have discretized space, our diffusion differential equation could be simplified as follows:
The table below provides an explanation for each variable.
Supposedly, the result of the fluorescent protein is supposed to be at least similar to a Turing pattern. However, because of the natural complexity of the Turing process in combination with an overflow of parameters, we were not able to generate a distinct Turing pattern that we desired. Before we tuned parameters and replacing terms in equations, our Turing patterns were not far from our diffusion result. The graph below shows our first run with the Turing patterning. There are some differences between this Turing plot and the diffusion result plot, but it is far from being significant enough for us to claim this modeling a successful replication of Turing patterning.
16 hours
We discussed with our math professors and decided to replace the terms that describes the degradation activators with the current term that we have. The modified result has improved in that it does appear more random. However, we decide that there is still no “pattern” within this ending result.
We ran our improved version of the model for 5, 10 and 16 hours. The comparison of the fluorescent protein concentration in those three time periods are shown below.
The result of this model did not achieve the Turing pattern that we have envisioned, but it has also demonstrated to us how difficult it is to form a completely stochastic Turing pattern. We will continue exploring this model either by modifying terms in the ordinary differential equation systems or by tuning the production and degradation rate of chemical substances. Although this is not a completely successful model, we learned how to improve our equations by looking at the biological circuit from a different perspective. We also learned that the slightest difference in parameters might cause a dramatic difference in the formation result, which leads us to appreciate the intricate nature of Turing patterning and appreciate its natural occurrence in nature.
Reference
Karig, D., Martini, K. M., Lu, T., DeLateur, N. A., Goldenfeld, N., & Weiss, R. (2018, June 26). Stochastic Turing patterns in a synthetic bacterial population.
Matlab code
Biofilm Lithography Model 1: Incorrect Hypothesis
Biofilm Lithography Model 2: Correct Hypothesis
Diffusion Model: chemical diffusion
Diffusion Model with Turing equations incorporated: Turing equations
Parameter values
All the numerical value of parameters that we have used: Parameters