A molecular dynamics simulation approach to model biofilm growth
Abstract
We developed a computational model to simulate the growth of biofilms on surfaces. The model applies methods from molecular dynamics and takes into account different physical and biological effects. In particular, interactions between bacteria and substrate- bacteria interactions are calculated with a force potential. The physical properties of the bacteria are updated using the Euler method. The biological properties are set using researched literature values as expected values for a normal distribution. To speed up the computation, we integrated multiprocessing. We evaluated the simulation results by comparing our results with recent literature. The model includes different sets of constants to model different bacterial strains. Furthermore, the software includes functions for visualization of the models' behavior over time. By providing an example notebook, well-documented code, and a python environment, we make the software easy to adapt for other teams' purposes. The model is published as a downloadable python package called "BiofilmSimulation".
Why did we do this?
The iGEM project of our team is the development of an "InToSens" (Inflammatory Toxin Sensor) made of biological components, which helps to diagnose the adherence of toxic biofilms to implant surfaces (click here for more Information). By detecting the attachment of a biofilm at an early stage, we aim to maximize the chance for treatment success. Therefore, our team is greatly interested in the mechanisms and parameters that cause and influence the formation of biofilms.
Because of the current pandemic, we are not able to research the biofilm in a wet lab. Hence, we decided to build a computational model to research the mechanisms of biofilm formation on implant surfaces.
The developed model takes into account a diverse set of biological and physical parameters and enables us, to quantitatively investigate the influence of different parameters on early biofilm formation.
We teamed up with the iGEM Team TU-Darmstadt to discuss our model throughout the development. The project of the Darmstadt Team aims to use biofilms to detox wastewater. Because of the lack of laboratory time, the Team is also interested in computational modelling of biofilm growth over time.
As both teams are interested in different biofilm compositions and influencing parameters, we made a great effort to make the model as flexible as possible.
In the following, we focus on the overview of the used methods and the models' capability. On our iGEM Git repository, we provide detailed instructions on how to use the model on your local machine. On this page, we take you through some of the considerations we made when building and implementing the numerical model.
How did we do this?
Throughout the developing process, we followed the principle of synthetic biology to design an artificial biological system by reducing the system on its main components.
We decided for OOP, because it provides a clear overview of the relationships between the different objects and the used parameters. In addition, an OOP implementation is easier to extend, for example by adding new interaction forces or even features in form of class methods.
This approach distinguishes our model from already existing ones, which often solve continuous equations or use cellular automatons ([1], [2]).
In detail, we provide the implementation of three classes and a variety of functions for visualization of the models' properties.
Used methods and implementation
Basic molecular dynamic scheme
For the simulation, we follow the basic scheme of MD simulations as depicted in Figure 1.
We initialize the simulation by setting initial parameters for each bacterium.
Subsequently, the program iterates over the initial bacteria and calculates the total force
on each bacterium by summing all of the acting forces. The index runs from 1 to N, the number of total bacteria in the biofilm.
\[ \vec{F}_i = \sum_i^N \vec{F}. \]
Based on the acting forces, the acceleration is calculated with the second Newton law of motion.
\[ \vec{a}_i = \vec{F}_i / m_i \]
The new position of the bacteria is calculated by using the Euler method.
Integrating the equation of motion results in the following relationship between velocity, acceleration and position \(r_i\)
\[ \vec{r}_i = \vec{v}_i \cdot \delta t + \frac{1}{2} \cdot \vec{a}_i \cdot {\delta t}^2 , \]
whereas the time step is represented in the formulas as \( \delta t \).
Afterwards the bacteria internal parameters e.g. size and mass are updated. After increasing the time step the process is repeated. We followed this scheme for our implementation.
Implementation
We provide the implementation of three classes and a variety of functions to visualize the simulation results.
Setting the initial conditions with the Constants class
The initial configuration of the biofilm can be specified by the user with an object of the Constants class. By creating an instance of the Constants class and changing its properties, simulation parameters like the simulation duration of the time step can be set. Furthermore, the number of initial bacteria and the bacteria strain can be set here and passed as an argument to the Biofilm. At the moment, two bacterial strains with the corresponding constants are available. The user can select between a biofilm consisting of Escherichia coli or Bacillus subtilis. Depending on the selected strain, different constants for the bacteria critical splitting length, mass, and the doubling time (and more) are used in the simulation. All of the used constants were researched in the literature and are listed in the functions' documentation and at the end of the page. The Constants class can be easily extended by other strains by adding the respective constants in the form of a python dictionary in the constants class.
Starting the simulation
After setting the simulation constants, the instance of the constants class can be passed to the Biofilm class, in which also the simulation is implemented. The Biofilm class is built from instances of the Bacteria class.
After setting the simulation constants, the instance of the constants class can be passed to the Biofilm class,
in which also the simulation is implemented. The Biofilm class is built from instances of the Bacteria class.
The simulation is started by spawning an initial configuration of bacteria on a plane.
These bacteria will be the origin of the biofilm in the simulation. All of the bacteria in the biofilm are stored in a
list as the parameter `bacteria_list`. By iterating over this list and updating the parameters of the bacteria objects inside the list multiple times,
the above MD simulation scheme (Figure 1) is realized. At the end of each iteration step, the current parameters of each bacteria object are saved and
appended to the parameters of the previous time steps.
Interaction Forces
The first calculations in each iteration steps are regarding the acting forces on each bacterium. We included the influence of gravitational force, the frictional force with the medium as well as Bacterium - substrate adhesion forces and Bacterium- Bacterium interaction forces.
The gravitational force is calculated by reading in the current mass of the bacteria and multiplying it with the gravitational acceleration on earth.
\[ F_g = m_\text{bacterium} \cdot g. \]
The frictional force with the surrounding medium is evaluated with the Stokes' law with the approximation for very small Reynold numbers and treating the bacteria as a spherical object.
The drag force is then calculated with
\[ F_\text{drag} = 6 \cdot \pi \cdot \mu \cdot R \cdot v \]
where \(F_d \) is the frictional force, \( \mu \) the dynamic viscosity of the medium, \( R \) the radius of the bacteria and \( v \) the bacteria drift speed.
In our simulation we use the bacteria length as the radius and the viscosity of the EPS matrix.
Furthermore, we decided to include two interaction forces in our simulation. The interaction of the bacteria with the substrate and surrounding bacteria is conducted by the Van-der-Waals force. Van-der-Waals forces occur for macroscopic bodies due to electronic dipole-dipole interactions of the molecules, of which the bodies are made of. Depending on the studied problem, different force potentials are used to model electronic interactions between macroscopic bodies.
We decided to use a Lennard-Jones potential of the form:
\[ V_\text{LJ}(r) = 4 \epsilon \left[ \left(\frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r}\right)^6 \right]. \]
Because we are interested in the resulting forces \( F_\text{LJ} \), we need to calculate and implement the gradient of the potential
\[ F_\text{LJ}( r ) = - \nabla V_\text{LJ} \]
which results in
\[ F_\text{LJ}(r ) = 48\epsilon \cdot \frac{\sigma^{12}}{r^{13}} - 24 \epsilon \cdot \frac{\sigma^6}{r^7} \]
This formula can be re-parametrized in terms \(r_\text{min}\), the distance at which the equilibrium is reached and the resulting force is \(0 \) and \(F_\text{min}\) representing the value of the global minimum.
This is accomplished by replacing \(\epsilon\) and \(\sigma\) with
\[ \sigma = \frac{r_\text{min}}{\sqrt[6]{2}} \quad \text{and} \quad \epsilon = F_\text{min}\frac{169\cdot r_{min}}{504 \cdot \sqrt[6]{7/13}}. \]
Implementing the re-parametrized in the function lennard_jones_force(r, f_min, r_min) is
a more convenient way than using \(\epsilon\) and \(\sigma\) as parameters for the force function, which can easily cause wrong values.
For the Bacterium- Bacterium interaction force, we set
\[ F_\text{min} = 6.81\; \text{nN}, \]
which is the maximal adherence force between bacteria measured in [3] and \(r_\text{min} = 2\) µm.
These values are use if the distance between the bacterias centre of gravity is greater more than \( 2.5 \) microns.
If the distance is smaller, we cut-off the Lennard-Jones potential and return the maximal force of \( 6.81 \) nN.
We do this to prevent overflow errors, which could occur when dividing by extremely small numbers.
For the Bacterium- Substrate interaction force, we set \(F_\text{min}= 7\) nN [3] and \(r_\text{min} = 1\) µm if the distance to the surface is greater than 5 microns. If the distance is smaller, \( F_\text{min} \) is returned. All of the forces are summed up and passed saved in the force parameter of each bacterium.
Updating the movement
We apply the Euler method to calculate the movement of the bacteria in three dimensions.
Using Newtons' second law of motion, the acting force is converted into an acceleration.
The velocity of the bacteria is then updated as
\[ \vec{v} = \vec{a} \cdot \delta t \]
and the position is calculated by integrating two times
\[ \frac{dv}{dt} = a \]
resulting in the new position
\[ \vec{r}_i = \vec{v}_i \cdot dt + \frac{1}{2} \cdot \vec{a}_i \cdot \delta t^2 \]
with the time step \( \delta t\).
Additionally, we add random movement to the bacteria by drawing values from a normal distribution with the expected value at the current position. By this, we include the influence of Brownian motion, which plays an important role in the early stages of biofilm formation [2].
Splitting and growing
Here we use the researched biological parameters to simulate the growth of the bacterium. The critical length serves as the expected value for a normal distribution, which decides if the bacterium is long enough to split. Furthermore the bacterium is grown according to the doubling time of the bacteria strain.
After these steps, the parameters are exported and stored in the JSON format before the next iteration step begins.
Constants used in the numerical model
The simulation is based on only a few constants, which we all researched in the literature and documented in our soaurce code. By collecting all of the used constants in one Class, we provide clear overview of which constant is used for which calculation. We guarantee, that no other constants are used except for the values listed in table 1-3.
We note, that some of the values serve as expected values for a normal distribution from which the actual value is drawn. The variances for the distributions are also researched and listed in the table.
| Parameter | Unit | expected value | variance | reference |
|---|---|---|---|---|
| Maximum bacteria substrate adhesion force | N | 7 * 1E-9 | 0 | [1] |
| Maximum bacterium- bacterium adhesion force | N | 6.81 * 1E-9 | 0 | [1] |
| Effective viscosity of EPS | Pa * s | log(1E5) | 0 | [2] |
| Effective viscosity of water | Pa * s | 0.7805 * 1E-3 | 0 |
| Parameter | Unit | expected value | variance | reference |
|---|---|---|---|---|
| Doubling time | s | 4200 | ||
| Growth rate | um / s | 5.239 * 1E-4 | [4] | |
| Division length | um | 4.7 | 0,564 | [4] |
| Birth length | um | 2.5 | 1 | [4] |
| Width | um | 1 | 0 | [4] |
| Mass | kg | 1,00E-12 | 0 |
| Parameter | Unit | expected value | variance | reference |
|---|---|---|---|---|
| Doubling time | s | 1200 | [6] | |
| Growth rate | um / s | 8.3 * 1E-4 | [6] | |
| Division length | um | 3.5 | 12,00 % | [5] |
| Birth length | um | 1 | 0.14 | [5] |
| Width | um | 0.5 | [4],[5] | |
| Free mean speed | um/s | 8 | [7] | |
| Mass | kg | 1,00E-12 | [4] |
Multiprocessing
With increasing bacteria numbers the models' computational cost increases very quickly.
Especially the calculation of interaction forces between the bacteria in the biofilm is lavish because
it is necessary to iterate over all bacteria two times. The order of complexity for the function bac_bac_interaction
is therefore \( O(n^2) \). Thus the usage of multiprocessing to increase the computational speed is inevitably necessary.
We make use of the multiprocessing package, which is a standard python package. For this, we rewrote all of the functions,
which change the parameters of the bacteria instances, to be able to apply the Pool class of the multiprocessing package.
With this class, we can map a function to the bacteria list and distribute the calculation on many processes,
resulting in a parallelized execution of the mapped function.
Results
As part of our Collaboration with the TU Darmstadt iGEM Team (see more here) we were provided SSH access to a Linux server with a CPU consisting of 32 cores. This enabled us to fine-tune the parameters with reasonable computation times. Furthermore, the TU Darmstadt Team adopted our Biofilm model for their purpose and verified, that the model is working as intended in numerous simulation runs with the Bacillus subtilis strain (Link to TU Darmstadt).
We applied our model for the Escherichia Coli strain.
In the following, we present the results of an exemplary simulation run with a duration of 360 minutes with time steps of 120 seconds. Because we were not able to grow and investigate a biofilm in the lab, we had to rely on experience from the literature to evaluate the results of our model.
The initial bacteria were placed randomly inside the range of 200 - 400 microns for the x and y-axis. The number 12 of bacteria increased significantly.
The plot in Figure 4 shows the number of living bacteria over the simulation duration. The number of bacteria increases exponentially with a generation time of 2391 seconds. This suits the doubling time of Escherichia Coli, which is reported in the literature to be around 1400 s long [6]. The end number of bacteria in the simulated biofilm is around 3000. This shows, that our implemented model can even process large bacteria numbers in a reasonable duration.
Despite we plugged in the constant "Doubling time",
the result is remarkable. This is because of the critical length, as well as the
birth length and mother-daughter size ratio are drawn independently from a normal distribution.
This can be verified with Figure 5, which shows the end distribution of the bacteria lengths in the biofilm.
The distribution expected value is around 2.2 microns. This fits the distributions for E. Coli. Length in the literature
[6].
Our simulation, therefore, models the bacterial growth pretty accurately.
The plot in Figure 5 shows the end positions of the bacteria on the surface.
One can see, that the bacteria moved away from their original position during the simulation.
The maximal movement range in 6 hours simulation time is about 0.8 mm. We notice that the bacteria
assemble as a cluster, while only a few bacteria lie apart. The formation of bacterial clusters
is also frequently reported in the literaturs
[4], [5].
Despite we simulated the emergence of only one of these clusters one can easily draw the conclusion, that on
medical implant surfaces many of these clusters will form rapidly. Only a few bacteria (in our case 12) are needed
to form an adherend bacterial cluster.
To confirm that modeling of the interaction forces works as intended, we started multiple simulations runs and switch on the interaction forces one by one. Each force contributes to the total force depicted in Figure 6. While at the beginning of the simulation, the forces on each bacteria do not vary much, we can see clearly that this changes with increasing simulation time (Figure 6). This is because in the early phase with small bacteria numbers, adhesion forces with the surface and the gravitational force prevail. With increasing bacteria numbers, the bacterium- bacterium interaction forces become more relevant. These interaction forces tend not to be fixed and cause many different values for the acting force.
The interaction forces range from very small numbers to 20 nN. We can report, that many bacteria with low force values probably correspond to the bacteria reaching an equilibrium in the interaction with the surface. The main contribution to the mean force, at about 7 nN, arises from the interactions between bacteria.
This range fits to the values reported in [7] and [8], which classify the adhesion forces on biomaterial surfaces from weak adhesion ranging from 0 to 2.5 nN, over strong adhesion (2.5 - 10 nN) to deactivated adhesion (> 10 nN). We were able to replicate these values with the numerical by using force potentials to model the interactions. The only fixed values we implemented, are regarding the form of the potentials. This is an awesome result and shows that the biophysical approach of our model is indeed working and modeling reality.
The bacteria velocities distribution approximates a normal distributed and reach values up to 0.6 µm / s (Figure 8). We notice, that the distribution is not entirely symmetric and the overall velocity is comparable small. The small absolute values of the velocities are due to the adhesion to the surface. We assume that with increasing simulation time more bacteria will show higher velocities because the bacteria distance to the surface will increase. The asymmetry is due to the drag force, which increases at higher velocities causing the bacteria to slow down and contributing to the distribution at lower values. We also observe a small peak at 0.1 µm / s. This could correspond to immobilized bacteria in the center of the biofilm.
Our simulation results are also shown as a 2D animation in Figure 9.
In this simulation, you can observe the behavior of the biofilm growth over 7 hours. The points
correspond to the position of the center of mass for each bacterium over time.
In the video, all of the above aspects can be observed.
We can see that the bacteria exhibit Brownian motion and spread on the surface.
We can see the formation of a cluster and watch the
interaction of bacteria with each other. For us the animation is a great way,
to review our numerical model. In combination with the plots above,
one can get a great impression of how the model behaves over time.
In conclusion, over biophysical molecular dynamics simulation models the behavior of biofilm growth
in the early phase astoundingly good. Despite rigorously reducing the very complex mechanisms of biofilm growth on a
few interactions and biological influences, we obtained great results, which are in line with the literature. We could confirm all of
the functionalities of our model by extensive testing and the iGEM Team Darmstadt even adopted our model for their purpose and carried out a
simulation with the Bacillus subtilis strain ( TU Darmstadt Biofilm Model).
Bibliography and references
- Yamamoto, Takehiro and Ueda, Shuya. 'Numerical Simulation of Biofilm Growth in Flow Channels Using a Cellular Automaton Approach Coupled with a Macro Flow Computation'. 1 Jan. 2013 : 203 - 216.
- Qin, B., Fei, C., Bridges, A. A., Mashruwala, A. A., Stone, H. A., Wingreen, N. S., & Bassler, B. L. (2020). Cell position fates and collective fountain flow in bacterial biofilms revealed by light-sheet microscopy. Science.
- Fang, H. H., Chan, K. Y., & Xu, L. C. (2000). Quantification of bacterial adhesion forces using atomic force microscopy (AFM). Journal of microbiological methods, 40(1), 89-97.
- JAMAL, Muhsin, et al. Bacterial biofilm and associated infections. Journal of the Chinese Medical Association, 2018, 81. Jg., Nr. 1, S. 7-11.
- RÖMLING, U.; BALSALOBRE, C. Biofilm infections, their resilience to therapy and innovative treatment strategies. Journal of internal medicine, 2012, 272. Jg., Nr. 6, S. 541-561.
- VAN HEERDEN, Johan H., et al. Statistics and simulation of growth of single bacterial cells: illustrations with B. subtilis and E. coli. Scientific reports, 2017, 7. Jg., Nr. 1, S. 1-11.
- ALAM, Fahad; KUMAR, Shanmugam; VARADARAJAN, Kartik M. Quantification of adhesion force of bacteria on the surface of biomaterials: techniques and assays. ACS Biomaterials Science & Engineering, 2019, 5. Jg., Nr. 5, S. 2093-2110.
- BUSSCHER, Henk J.; VAN DER MEI, Henny C. How do bacteria know they are on a surface and regulate their response to an adhering state?. Plos pathog, 2012, 8. Jg., Nr. 1, S. e1002440.