Team:IISER Kolkata/Biofilm Degradation
![[preloader image]](https://static.igem.org/mediawiki/2021/d/d6/T--IISER_Kolkata--image-preloader.gif)
Loading . . .
Bacteria form biofilms as part of their survival mechanisms, and biofilms are thus ubiquitous. The biofilms are clusters of bacteria that are attached to a surface and/or to one another and embedded in a self-produced matrix. The biofilm matrix consists of substances like proteins (e.g., fibrin), polysaccharides (e.g., alginate), as well as eDNA. In addition to the protection offered by the matrix, bacteria in biofilms can employ several survival strategies to evade the host defence systems.
In this mathematical model, we try to understand how the biofilm formed by the pathogen in an infected cow reacts in the presence of the chosen biocide in our project, Nisin. We try to visualise how this biofilm dome-like structure ruptures and falls apart on treatment with Nisin in a 2-D capillary tube.
Here we consider a 1mm \(\times\) 1mm two dimensional grid segment of the capillary tube horizontally centered around the biofilm. The biofilm is a semicircle attached to the bottom of this grid setup with its centre at 0.5mm from the left end and its radius 0.3mm. This biofilm consists of four layers of with a varying degree of live and dead bacteria (initially there is no dead bacteria in these layers). This biofilm dome is treated with a flowing solvent containing the antimicrobial agent of choice (here it is Nisin). The fluid dynamics follows Navier- Stokes equation with the added condition of incompressiblity. We map this segement of flow into a 128 $\times$ 128 matrix-grid with all the necessary parameter defined is a their corresponding matrices of the same dimensions.
- The growth of the biofilm in the treatment phase is negligible.
- The extra cellular polymeric substances are absent.Film only consists of live and dead bacteria.
- Biocide penetrates the biofilm mainly by diffusion.
The reaction–diffusion–advection equations for d, \(\phi_{bl}\), \(\phi_{bd}\) :
\begin{array}{l} \frac{\partial}{\partial t}\left(\phi_{s} d\right)+\nabla \cdot\left(\phi_{s} \mathbf{v} d-D_{e} \phi_{s} \nabla d\right)=-\phi_{b l} \frac{K_{1} d}{d_{0}+d}\rightarrow (1) \\ \frac{\partial \phi_{b l}}{\partial t}+\nabla \cdot\left(\phi_{b l} \mathbf{v}\right)=\nabla \cdot\left(\Lambda \phi_{b l} \nabla \frac{\delta f}{\delta \phi_{b l}}\right)-\phi_{b l} \frac{K_{2} d}{K_{d}+d}\rightarrow (2) \\ \frac{\partial \phi_{b d}}{\partial t}+\nabla \cdot\left(\phi_{b d} \mathbf{v}\right)=\nabla \cdot\left(\Lambda \phi_{b d} \nabla \frac{\delta f}{\delta \phi_{b d}}\right)+\phi_{b l} \frac{K_{2} d}{K_{d}+d}\rightarrow (3) \end{array}Chemical free energy density equation:
$$ f\left(\phi_{b}\right)=k_{B} T\left[\frac{\Gamma_{1}}{2}\left\|\nabla \phi_{b}\right\|^{2}+\Gamma_{2}\left(\frac{1}{N} \phi_{b} \ln \phi_{b}+\left(1-\phi_{b}\right) \ln \left(1-\phi_{b}\right)+\chi \phi_{b}\left(1-\phi_{b}\right)\right)\right]\rightarrow(4) $$Chemical potenital equation:
$$ \frac{\delta f}{\delta \phi_{b}}=k_{B} T\left[-\Gamma_{1} \Delta \phi_{b}+\Gamma_{2}\left(\frac{1}{N}\left(1+\ln \phi_{b}\right)-\ln \left(1-\phi_{b}\right)-1+\chi\left(1-2 \phi_{b}\right)\right)\right] $$from (4) and by adding (2) and (3) we get,
$$ \frac{\partial \phi_{b}}{\partial t}+\nabla \cdot\left(\phi_{b} \mathbf{v}\right)=\nabla \cdot\left(\Lambda \phi_{b} \nabla \frac{\delta f}{\delta \phi_{b}}\right) $$Modified Navier Stokes Equation:
$$ \begin{array}{l} \nabla \cdot \mathbf{v}=0 \\ \rho\left(\frac{\partial \mathbf{v}}{\partial t}+\mathbf{v} \cdot \nabla \mathbf{v}\right)=\nabla \cdot(v \mathbf{D})-\nabla p+\frac{\delta f}{\delta \phi_{b}} \nabla \phi_{b} \end{array} $$with:
$$ v=v\left(\phi_{b d}, \phi_{b l}, \phi_{s}, v_{b}, v_{s}, d\right), \quad \rho=\phi_{b} \rho_{b}+\phi_{s} \rho_{s}, \quad \mathbf{D}=\frac{1}{2}\left(\nabla \mathbf{v}+\nabla \mathbf{v}^{T}\right) $$Equation for Effective Diffusion Coefficient:
$$ D_{e}=\frac{2 D_{d}^{0}\left(1-0.4{\phi}_{b}\right)}{\left(2+0.4{\phi}_{b}\right)\left(\phi_{s}+0.6{\phi}_{b} / D_{p r}\right)} $$Effect of viscosity of the material:
$$ v=\phi_{b} v_{b}+\phi_{s} v_{s} $$Numerical methods:
$$ \tilde{\mathbf{x}}=\frac{\mathbf{x}}{h_{0}}, \quad \tilde{t}=\frac{t}{t_{0}}, \quad \tilde{\mathbf{v}}=\frac{\mathbf{v} t_{0}}{h_{0}}, \quad \tilde{p}=\frac{p t_{0}^{2}}{\rho_{0} h_{0}^{2}}, \quad \tilde{d}=\frac{d}{d_{0}} $$| Symbol | Parameter | Value | Unit |
|---|---|---|---|
| \(T\) | Temperature | 296 | Kelvin |
| \(\Gamma_{1}\) | Distortional energy | \(8 \times 10^{6}\) | \(\mathrm{~kg} \mathrm{~m} \mathrm{~s}^{-2}\) |
| \(\Gamma_{2}\) | Strength of the bulk mixing energy | \(3 \times 10^{17}\) | \(\mathrm{~kg} \mathrm{~m} \mathrm{~s}^{-2}\) |
| \(\chi\) | Flory-Huggins mixing parameter | \(0.55\) | dimensionless |
| \(N\) | Generalized polymerization parameter | \(1 \times 10^{3}\) | dimensionless |
| \(\Lambda\) | Mobility parameter | \(1 \times 10^{-9}\) | \(\mathrm{~kg}^{-1} \mathrm{~m}^{3} \mathrm{~s}\) |
| \(v_{b}\) | Viscosity of the bacteria | 5 | \(\mathrm{~kg} \mathrm{~m}^{-1} \mathrm{~s}^{-1}\) |
| \(v_{s}\) | Viscosity of the solvent | \(1.002 \times 10^{-3}\) | \(\mathrm{~kg} \mathrm{~m}^{-1} \mathrm{~s}^{-1}\) |
| \(\rho_{b}\) | Biofilm density | \(1 \times 10^{3}\) | \(\mathrm{~kg} \mathrm{~m}^{-3}\) |
| \(\rho_{s}\) | Solvent density | \(1 \times 10^{3}\) | \(\mathrm{~kg} \mathrm{~m}^{-3}\) |
| \(D_{N}^{0}\) | Aqueous diffusion coefficient of Nisin | \(0.19 \times 10^{-9}\) | \(\mathrm{~m}^{2} \mathrm{~s}^{-1}\) |
| \(K_{1}\) | Maximum biocide consumption rate for nisin | \(2 \times 10^{-7}\) | \(\mathrm{~kg} \mathrm{~m}^{-3} \mathrm{~s}^{-1}\) |
| \(d_{0}\) | Biocide Monod half-saturation coefficient | \(1 \times 10^{-3}\) | \(\mathrm{~kg} \mathrm{~m}^{-3}\) |
| \(K_{2}\) | Maximum decay rate of live bacteria (due to biocide action) | \(2 \times 10^{-5}\) | \(\mathrm{~s}^{-1}\) |
| \(K_{d}\) | Live bacteria decay Monod half-saturation coefficient | \(1 \times 10^{-5}\) | \(\mathrm{~kg} \mathrm{~m}^{-3}\) |
| \(\bar{d}\) | Reference biocide concentration for weakening biomaterial | \(5 \times 10^{-3}\) | \(\mathrm{~kg} \mathrm{~m}^{-3}\) |
| \(\bar{T}\) | Reference time scale for weakening biomaterial | 120 | \(\mathrm{~s}\) |
| \(h_{0}\) | Characteristic length scale | \(10^{-3}\) | \(\mathrm{~m}\) |
| \(t_{0}\) | Characteristic time scale | \(5 \times 10^{-2}\) | \(\mathrm{~s}\) |
| \(\rho_{0}\) | Characteristic density scale | \(1 \times 10^{3}\) | \(\mathrm{~kg} \mathrm{~m}^{-3}\) |
The following parameter matrices need to be initialised before running the simulation:
- Volume fraction of live bacteria
- Volume fraction of dead bacteria
- Volume fraction of total bacteria
- Volume fraction of solvent
- Concentration of Nisin
- Horizontal and Vertical components of fluid velocity in the system
- Density across the system
- Effective diffusion Coefficient
- Pressure across the system
- Viscosity across the system
- Chemical potential
Setting up the Initial Velocity of the fluid
Horizontal component of velocity
The Horizontal component of velocity of solution at every point will be symmetric about the vertical line passing through the centre of the biofilm. Hence the velocity profile of the left half is replicated in the right half.
The first and the last of the 128 columns has a quadratic velocity profile with bottom-most layer having zero velocity and top layer having a maximum speed of 2cm/s.
The column through the centre of the biofilm has zero horizontal velocity in the innermost layer of biofilm and maximum velocity at the top. The horizontal velocity variation in this column is bi-quadratic.
For the horizontal component of velocity in all the rows above the biofilm From the equation of continuity we see that for a particular row, speed of the fluid increases as we move towards the centre. This is because height available for the fluid decreases (due to presence of biofilm) as we move towards the centre. So our velocity profile will have a maximum in the middle and minimum at the ends for a given row in the grid map. Because the incompressibility condition has to be satisfied, these three points have to be turning points.The simplest variation that has three turning points (leftmost, rightmost and at the middle) is a bi-quadratic. The horizontal component of velocity in all the rows above the biofilm will have bi-quadratic variation symmetric about the middle.
In the middle, right above the centre of the biofilm using equations of continuity and the fact that the velocity will be purely horizontal at the top of the dome of the biofilm .There will be a quadratic variation of this horizontal velocity from the top of the solvent flow to the innermost layer of the biofilm, with the innermost layer of biofilm having 0 velocity and top layer having maximum speed.With this, we calculate a new maximum speed at this column location.
In gridpoints to the left and right of the biofilm In a particular row below the height of 0.3mm from the bottom, the velocity varies in a bi-quadratic fashion from the leftmost gridpoint having the maximum velocity to the point in the innermost interface of biofilm having zero velocity. This can be extended to right of the biofilm due to symmetry.
Vertical component of velocity
The vertical component of velocity of solution at every point will be symmetric in magnitude about the vertical line passing through the centre of the biofilm but opposite in direction i.e, in the left half the velocity direction is upwards while in the right half its downwards. Hence the velocity profile of the left half is replicated in the right half in magnitude and reversed in direction.
The vertical velocity of the innermost layer of the biofilm is obtained by multiplying the slope of inner layer at each grid point to the horizontal velocity component.
The vertical velocity of the other outer layers of the biofilm is calculated from the incompressibility of the fluid condition.
In other regions the vertical velocity varies quadratically.
Setting up concentration grid
Bacterial and solvent concentration: The biofilm has four layers. The outermost layer has bacterial volume fraction of 0.1, the layer beneath has 0.2,the layer beneath the second 0.3 and the innermost layer has 0.4. Inside these four layers the bacterial volume fraction is 1. Everywhere else in the grid the bacterial volume fraction is zero.
Solvent volume fraction= 1-bacterial volume fraction.
Initially there is no dead bacteria, and bacterial volume fraction is live bacterial volume fraction.
Antimicrobial Concentration(matrix 'd'): Same throughout the grid
Matrix initialisation for viscosity, density, pressure and effective diffusion coefficient can be easily done using governing equations.
We first simulated the differential equations using the parameters reported in the literature (as given in the parameters table). This involved simulating the system in the presence of a biocide with known parameters so that we feel confident about the reproducibility of the model by extension, the validity of the results. We then simulated the systems using parameters obtained from molecular dynamics simulations.
Spatial behaviour
The above figure shows the progression of the volume fraction of dead and live bacteria with simulation time. The first row is the initial volume fraction of bacteria (at T=0). The second row is the volume fraction of live bacteria at each T. The third column represents the volume fraction of dead bacteria at each time T. The X and Y axes represent grid cells. We have divided a 1 sq. mm surface into 128x128 spatial unit cells.
Here is an animation to better elucidate the behaviour of the biofilm that ultimately leads to the subsequent degradation of biofilm.
It is clear from the animation that the biofilm degrades completely after a stipulated time. The degradation of the biofilm begins at the top of the biofilm dome and progresses down the Y-axis. This is natural considering that the outermost layer of biofilm is the first line of defence against biocide.
Comparing the efficacy of NisinA and NisinPV
This is the part of the model that links the biofilm degradation model at the tissue level to the molecular level as well as the population-level model.
From the molecular dynamics simulations, we got quantitative values for the diffusion coefficient of both NisinA and NisinPV. Additionally, from the free energy of binding to the NSR molecule, we were able to make an estimate of the kinetics of binding using the relationship between free energy and equilibrium constant:
$$\Delta G^{\circ}=-R T \ln K $$ $$K=e^{-\frac{\Delta G^{\circ}}{R T}}$$Considering that the free energy of binding to NSR for NisinA is 0.8775 kcal/mol and for NisinPV is 4.6292 kcal/mol, at equilibrium, the rate of binding of Nisin PV will be 6.4 times lower than that of Nisin A (at 295K). The binding rate is proportional to the rate of degradation of nisin. Therefore, we were able to distinguish between the relative rate of degradation of Nisin A and Nisin PV. We simulated the differential equation set up in the biofilm degradation model with these values of the diffusion coefficient and degradation constants. We simulated the systems over an array of concentrations of each biocide and plotted the average bacterial volume fraction across the model space with time.
The results of this simulation are outlined in the following animations.
From these animations, it can be inferred that under the given conditions, the average bacterial dead bacterial volume fraction reaches a saturation point at a significantly lower concentration of biocide in the case of NisinPV as compared to NisinA (approximately two orders of magnitude difference).
| Biocide concentration (mM) at which saturation is observed | Corresponding maximum dead bacteria volume fraction | |
|---|---|---|
| Nisin A | ~ 1.5 E(-4) | ~ 0.3.5 |
| Nisin PV | ~ 1.9 E(-6) | ~ 0.45 |
Given the complexity of the simulation and computational limitations, we could not simulate the system for longer than 1000 time units and therefore we could not explore the behaviour at larger timescales. Had we been able to simulate the system for longer, we would have quickly arrived at the concentration at which the average dead bacteria volume fraction would have approached 0.5. We could consider the concentration corresponding to the volume fraction of dead bacteria reaching 0.5 as a measure of the inhibitory concentration of the biocide. The IC-50 value of Nisin is reported in the literature in the nano-molar range. We can see that Nisin PV concentrations are already in the nanomolar range.
Having an idea about the inhibitory concentration range of Nisin PV, we can implement the same in terms of population dynamics and consequently human practices.
Assuming that the trend of bacterial volume fraction continues, it can be assumed that the IC-50 will be reached at around 10 nM. From the gene circuit model, we calculated the rate of production of Nisin to be in the order of 10^(-4) uM per bacteria. Therefore a typical inoculum size of a few million bacteria per mL should suffice to reach the IC-50 value within 20 hours. Also, we have DNAse-I in the system in our actual GMO, which would accelerate the biofilm degradation. Therefore, we are making conservative estimates.
References
- Zhang T, "Modeling of biocide action against biofilm," Bulletin of Mathematical Biology, Montana State University.