Individual based Modelling of Bioﬁlm Colonisation

: Bacterial bioﬁlms on implants or host tissue can cause unwanted infections. One of the main factors for the pathogenic bioﬁlm growth is the initial adhesion of the bacteria. In this project a framework for systematic studies of the bacterial adhesion is provided. Therefore, an Individual based Model (IbM) is set up in the IbM solver NUFEB, an extension of the molecular dynamics code LAMMPS, which integrates Newton’s equations of motion for the cells that interact via various forces. Experimentally measured adhesion forces are implemented for the cell-wall and cell-cell interaction. Furthermore, the random movement of planktonic cells is included into the IbM software. A sensitivity analysis of important parameters and an exemplary simulation of the attachment are performed. Thus, the attachment of the bacterial cells and the initial growth of a bioﬁlm can be qualitatively modelled.


Introduction
Biofilms are structured bacterial communities which are attached to a surface and are held together by selfproduced extracellular polymeric substances (EPS) [8]. They are the prevailing habitat of bacteria and protect the individual microorganism from starvation, desicca-tion, the host immune system or drug treatment [6,7]. The formation of a biofilm consists of several steps and is influenced by many factors such as temperature, pH value, osmotic interactions, mixing and oxygen concentration [1]. First, floating planktonic bacteria have to be transported to the surface. They have to overcome e.g. certain surface conditions like surface charge or roughness, fluid flow, Coulomb forces or hydrophobic interactions to attach to the surface [17,29]. This can be achieved via Brownian motion, sedimentation, convective mass transport by the movement of the bulk fluid or by active movement of the bacteria [29]. Often smaller microcolonies consisting of several layers are formed and the bacteria start to multiply [17]. Bacterial cells can communicate via quorum sensing by releasing small diffusible molecules, the so-called autoinducers [32]. When a certain density is reached a synchronised change in the gene expression takes place facilitating the secretion of EPS [26]. The biofilm reaches a maximum growth where cell division and death are in equilibrium [1]. Cells can be released by sloughing due to mechanical stress caused by the fluid flow, active dispersion or lost from predation and viruses [19,30]. Finally, the breakdown of the biofilm occurs, where bacteria detach and return into the planktonic state [1]. Biofilms can be beneficial in for example waste water treatment but can cause problems in the medical context where they cause infections in the host tissue or on implants and in technical applications where they induce biofouling or biocorrosion [5,13,23]. Bacterial biofilms cause 80% of human infections [4]. For example in the USA biofilm-associated diseases are responsible for over half a million deaths annually and cause immense annual treatment costs of about $94 billion [33]. Biofilms that induce diseases are termed pathogenic biofilms. Adhesion is the essential step in the pathogenesis with bacteria being killed by the host immune system if they do not attach quickly [18]. Pathogenic biofilms occur for example in the oral cavity where many different bacteria are found. First, initial colonisers e.g. streptococci (e.g. Streptococcus oralis) and actinomyces (e.g. Actinomyces naeslundii) adhere to the surface of teethes and implants. These pioneer colonisers are not pathogenic but provide attachment areas for secondary and late colonisers allowing pathogens e.g. Porphyromonas gingivalis to settle so that a multispecies biofilm forms [9]. The biofilm environment provides advantages for the individual microorganism compared to the planktonic lifestyle. It increases protection against antibiotic treatment, but also against extreme environmental conditions, phagocytosis and dehydration which in return increases the pathogenicity of the individual bacterium [7,11]. Due to the secretion of EPS the biofilm is more resistant against antibiotics and the amount of antibiotic needed can be 100 − 1000 times more than for planktonic bacteria [1]. In order to avoid infections, the formation, especially the adhesion of bacteria has to be prevented. For a long time, the only method to treat the biofilm formation was increasing antibiotics or by removing the infected tissue/device. Nowadays, novel materials such as nanomaterials for implants are used to prevent pathogenic biofilms [4]. In addition to experiments, mathematical modelling of biofilms offers a deeper insight into the process of biofilm formation to prevent infections and efficiently use nonpathogenic biofilms. Many mathematical models exist, varying in complexity in treating the biomass and its transport. Mainly two different approaches have been followed: continuum models, where the bacteria are treated as concentrations, and discrete models, where the bacteria are treated as individuals. The domain is usually divided into the bulk liquid, the boundary layer and the biofilm [25]. The models include a transport mechanism like diffusion and advection, a consumption and a growth mechanism as well as a loss mechanism [19]. Continuum biofilm models started to develop in the 1970s starting with one dimensional formulations and expanding to complex multispecies multidimensional dynamic models. In the 1990s discrete biofilm models started to arise developing the biomass representation and its spreading mechanisms. Three different groups have been classified: Cellular Automaton models (CA), hybrid differential-discrete CA models and Individual based Models (IbM) [25]. IbMs of biofilms have been introduced by Kreft et al. [21] and bear the advantage that the individual characteristics and interactions of the bacteria are considered and the emergent behaviour of the biofilm can be analysed. They can be used to analyse the adhesion process in detail. Different IbM codes for the modelling of biofilm growth have been introduced: BacSim [21], iDynoMiCS [22], Simbiotics [28] and NUFEB [24]. Koshy-Chenthittayil et al. [20] review the different models in their article. Adhesion is the essential factor for biofilm development which can cause infections in the presence of pathogenic bacteria when they attach to implant surfaces or to the host tissue. Nevertheless, most biofilm models neglect these early stages of biofilm formation and focus on the mature biofilm growth and detachment. Achinas et al. [1] describe the bacterial adhesion as one of the remaining challenges in microbial research. The goal of this project is to provide a foundation for systematic studies on influences on the adhesion process in order to develop new combat strategies. In this project, a biofilm colonisation model is set up using Individual based Modelling. Thus, the adhesion process can be described on bacterial level. Most of the available IbM solvers for simulation of the 3D dynamics of microbial communities focus on the biological aspects of biofilm growth but the used open-source IbM solver NUFEB 1 includes also mechanical and chemical processes and offers coupling with fluid dynamics [24]. NUFEB is an extension to the general purpose opensource C++ molecular dynamics code LAMMPS 2 . The remainder of the article is structured as follows: In Section 2 an overview of NUFEB is given and the extensions thereof are explained. The bacterial adhesion to the surface and the cell-cell adhesion are modelled via experimentally measured adhesion forces. Furthermore, the code is extended to include the random movement of planktonic cells. In Section 3 the results of several simulations are depicted to examine the attachment of bacterial cells. Finally, in Section 4 the results and the set-up of the IbM are discussed.

Force-distance curves
For the interaction between the bacteria and the wall measured adhesion forces have been provided by Dr. Katharina Doll (Hannover Medical School, Department of Prosthetic Dentistry and Biomedical Materials Science, Hannover, Germany). The force-distance curves have been measured using single-cell force spectroscopy with a FlexFPM atomic force microscope (Nanosurf AG, Liestal, Switzerland) connected to a FluidFM pressure control system (Cytosurge AG, Zürich, Switzerland) and mounted on an inverse microscope (Eclipse Ti-S, Nikon GmbH, Düsseldorf, Germany). A hollow silicon nitride cantilever with a circular opening of 300 nm at the pyramidal tip with a theoretical spring constant of 0.6 N/m was used. With the cantilever a single cell was targeted with a setpoint force of 10 nN and a negative pressure of 400 mbar for 5 s and retracted with a velocity of 1 µm/s. The bacterium was then transferred to a titanium surface with a setpoint force of 0.75 nN, an adhesion time of 5 s, a velocity of 1 µm/s with enabled force feedback. This measuring method is described in Doll et al. [9]. A schematic force-distance curve can be seen in Figure 1 from which the adhesion forces and detachment distance can be identified. The adhesion forces have been measured for different species e.g. S. oralis and P. gingivalis.

Implementation in LAMMPS
LAMMPS is a classical molecular dynamics code used to model particles in a liquid, solid or gaseous state by providing the possibility to include various interactions and boundary conditions. NUFEB extends the code for microbial simulations. In LAMMPS a "fix" defines an operation that applies during the simulation. Following a short overview of NUFEB they will be subsequently defined for the new operations.

Overview of NUFEB
A short overview of the IbM solver NUFEB is given following the description of Li et al. [24]. For the description the ODD protocol (Overview, Design concepts and Details) introduced by Grimm et al. [12] is followed.

Model overview
The IbM's purpose is to model the heterogeneity and interactions of the individual microbes to study and predict the emergent properties and behaviours of e.g. the biofilm.
The microbes are represented as spheres with individual attributes like position, density, velocity, force, inner and outer mass, inner and outer diameter, growth rate, yield etc. which change over time. In the domain, a box with dimensions L X ×L Y ×L Z is defined and the chemical properties like nutrient concentrations, pH and Gibbs free energies are continuous and can be calculated at each discretised grid point.

Design concepts
The following concepts apply: 1. Emergence. From the individual attributes and interactions the morphology and spatial distribution at the population level emerge. 2. Sensing. The individual growth depends on the local chemical properties and the motion on the local mechanical interactions and fluid flow. 3. Observation. The physical, biological and chemical states are stored at every time step. 4. Interaction. The interaction between individuals result from mechanical forces. 5. Stochasticity. The initial distribution of microbes, the size and distribution of daughter cells are considered as stochastic processes.

Details
NUFEB consists of several sub-modules. There is the biological sub-module, where the growth and decay of the microorganisms are calculated by ordinary differential equations: where m i is the biomass and µ i is the specific growth rate of the bacterium i which can either be specified by a Monod-based or by an energy-based process. The Monod-based process is based on the local nutrient concentrations. The decay rate is of first order. Three different functional groups, active heterotrophs (HET), nitrite oxidizing bacteria (NOB) and ammonia oxidizing bacteria (AOB), are implemented along side inactive EPS and dead cells. For example the growth rate for the GAMMAS 2022 aerobic growth of HET is calculated as follows using the Monod-based approach: where µ max is the maximum specific growth rate of the biomass group, S S is the local concentration of the substrate, S O 2 the local concentration of oxygen, K S,HET the affinity constant between the substrate and the biomass group HET and K O 2 ,HET between oxygen and HET, respectively.
Division and death are treated as instantaneously where division takes place if a certain threshold is reached. Then, the mass is divided into two daughter cells randomly assigned 40% − 60% of the mother cell. One cell takes the position of the mother cell and the other cell is randomly located at a distance d . Microbes are considered dead if they shrink below a threshold diameter and cannot perform biological processes anymore, but shrink linearly. EPS production can be implemented for active heterotrophs using the Monod-based growth model. First, EPS is secreted as a shell around the individual. If a threshold is reached, a separate EPS particle is randomly placed next to it. The physical sub-module contains the interaction between the individual cells as well as with the surrounding fluid and uses the discrete element method (DEM). Mechanical relaxation is performed to establish the mechanical equilibrium which could have been altered due to the growth and division of the cells and the resulting collisions or overlapping. Therefore, the Newtonian equations of motion are solved for the movement of every cell i : where v is the velocity of the cell and F is the sum of all effective forces. For solving the Velocity Verlet algorithm is used which relies on the neighbour lists. Neighbour lists contain the neighbours which are located within a specified radius for every cell. In NUFEB, so far the following forces have been implemented: F c is the contact force, F a is the EPS adhesive force and F d is the drag force which are all frequently used forces in microbial systems. The contact force is based on Hooke's law: where N i is the total number of neighbouring cells of the cell i , K n is the elastic constant for normal contact, δn i , j is the overlap distance between cells i and j , m i , j is the effective mass of cells i and j , γ n is the viscoelastic damping constant for normal contact and v i , j is the relative velocity of the two cells.
The EPS adhesive force is calculated as a Van der Waals force: where H a is the Hamaker constant, r i , j is the effective outer radius of cells i and j , h min,i , j is the minimum separation distance of the two cells and n i , j is the unit vector between i and j . Finally, the drag force is computed as: where f,i is the cell volume fraction with f,i = 1 − s,i , V p is the volume of the cell, u p its velocity, U f the fluid velocity acting on the cell and β the drag correction coefficient.
In the chemical sub-module the nutrient transport is described via the diffusion-advection-reaction equation. It is solved for each soluble concentration S: where D is the diffusion coefficient, U is the fluid velocity and R is the reaction term. The reaction rate is defined according to the Monod kinetics or the energybased model. It is solved using a discretised Marker-And-Cell grid where the concentration is defined at the center of a cubic. Temporal and spatial derivatives are discretised by forward Euler and central finite differences. Different boundary conditions as non-flux Neumann, Dirichlet or periodic can be chosen. NUFEB can be coupled with fluid dynamics using the CFD-DEM solver SediFoam, an interface between the solvers LAMMPS and OpenFOAM. Cell information are transferred to OpenFOAM where an averaging procedure is applied to convert them to the Eulerian CFD mesh where the fluid momentum equations are solved using the PISO algorithm. Then, the drag force and velocity field are transferred to LAMMPS again. The IbM simulation is summarised in Algorithm 1. First, if included, the fluid dynamics are solved and the resulting drag force and the new fluid velocity field are calculated (see Line 1.1). Next, the nutrient mass balance is solved to calculate the nutrient concentrations (see Line 1.2). The nutrient concentration in the bulk is updated (see Line 1.3). Then, the mass and size of the microbes are updated by performing the microbe growth (see Line 1.4), which is followed by their division, Algorithm 1 Pseudo code for IbM simulation in NUFEB (adapted from Li et al. [24]) 1: solve fluid dynamics and update velocity field and cell drag force 2: solve nutrient mass balance and update concentration field 3: update nutrient concentration in bulk based on total consumption rate in biofilm 4: perform microbe growth and update mass and size 5: perform microbe division, death and EPS production 6: perform mechanical relaxation and update locations and velocity 7: update boundary layer and neighbour list death and EPS production (see Line 1.5) and the performance of the mechanical relaxation where the different forces acting on the cell are considered to update the individual locations and velocities (see Line 1.6). Finally, the boundary layer and the neighbour lists are updated (see Line 1.7). For the coupling the processes are assumed as pseudo-steady state because the time steps are of different magnitudes (mechanical ≈1 × 10 −7 s, diffusion ≈1×10 −4 s and biological from minutes to hours).

Wall interaction
For the interaction between the wall and the bacteria, experimentally measured adhesion forces should be included when the cells have attached to the wall and are withdrawn by some force from the wall. In order to use the tabulated forces for this interaction, the code is extended. The syntax of the command can be seen In a next step, it is checked if the withdraw variable equals one. If true, the corresponding adhesion force is computed from the imported table via linear interpolation for the distance between the cell and the wall. The distance and corresponding force values from the table are shifted by the radius of the cell, so that the bacteria do not penetrate into the wall (see Line 2.8). Afterwards, it is checked whether the exerted force on the cell in the specified dimension which is perpendicular to the specified wall (x = 1, y = 2 and z = 3) is pointing away from the wall and larger than the tabulated force. If it is larger, the tabulated force is added to the overall force (see Line 2.10). If it is smaller, the overall force is set to zero so that the cell is not retracted from the wall and does not move (see Line 2.12). In the simulation the forces resulting from the cell-wall interaction are considered at last in the mechanical relaxation. Thus, all the relevant forces are considered in the simulation and the forces are set to zero so that the cell remains at the wall if the effective forces are not large enough to move the cell.

Cell-Cell interaction
For the cell-cell interaction a similar approach as for the cell-wall interaction is pursued.

This loop contains an inner loop over all neighbours
numberOfNeighbours of the cell currentCell. In the inner loop, for each neighbour currentNeighbour, the corresponding force from the imported table is computed by subtracting both cell radii from the cell center to cell center distance and calculating the corresponding force via linear interpolation (see Line 3.3). Next, the unit vector for the force n due to the adhesion is calculated pointing from the cell center of the current cell to the corresponding cell center of the neighbouring cell. Furthermore, the unit vector pointing in the direction of the overall force f is calculated. Then, the vector n is projected onto f to extract the fraction of the adhesion force that is pointing in the direction of the overall force f (see Line 3.4). The projection of n on f is calculated as follows: The corresponding forces are depicted in Figure 2. The resulting vector p is multiplied by the interpolated value from the table to obtain the correct effective adhesion force p currentNeighbour . Then, only this force is added to the total adhesion force p currentCell . Furthermore, it is checked if cells are overlapping. The overlap is checked, so that forces can be exerted, if overlaps between cells occur (see Line 3.7). The adhesion forces are only considered if no overlaps occur. After all forces from the neighbours have been added, in a next step, it is checked if the overall force f currentCell is larger than the total adhesion force p currentCell and no overlap occurs. The overall force is calculated from the cell-cell interaction to resolve the overlaps. Then, the adhesion force is added to the overall force (see Line 3.11).
If it is smaller and no overlap occurs, the overall force is set to zero and the cell does not move (see Line 3.13).
If an overlap occurs, the adhesion forces are not considered, so that first the overlaps are removed and in the next time step the adhesion forces can be taken into account. The cell-cell interaction are considered in the beginning of the simulation after the forces resulting from the Lennard-Jones potential have been considered. At later stages, external forces e.g. the drag force from the fluid flow are considered and can be added to the total force.

Planktonic behaviour
As a further extension of the code and to improve the modelling of the adhesion behaviour of bacteria the possibility of adding planktonic cells to the simulation is included. Therefore, a new fix planktonic is added, for which the command is displayed in Listing 3, where ID specifies the name of the fix and group-ID specifies the group of bacteria the fix should be applied on. to the biofilm-associated type and cannot perform a random movement to simulate the attachment of planktonic cells (see Line 4.14).

Individual based Model
To analyse the attachment behaviour of the bacteria an exemplary IbM simulation for the species S. oralis with the implemented operations is set up. The domain is initialised as a box with lower and upper limits (136 µm × 136 µm × 20 µm). The x and y boundary are periodic and the z boundary is fixed. For this simulation two cell groups are defined, the predefined heterotrophic (HET) cells and the group of planktonic cells (pla). The HET cells are the biofilm-associated cells that perform growth and division and the planktonic cells perform a random movement. The process of the IbM simulation is summarised in Algorithm 5. First, the parameters for the chemical processes are set, which are updated every 100 time steps as there are 5760 time steps for the exemplary simulation and can be updated accordingly as the conditions do not change that fast. The domain is divided into 68 grid elements in x and y dimension and 10 into the z dimension. The diffusion time step is set to 1 × 10 −4 s. The maximum number of iterations in the kinetics integration is set to 5000. The mass balance of the soluble substrates (see Line 5.1) is solved with a defined absolute tolerance for convergence of 1 × 10 −6 and the diffusion factor is defined as constant (see Equation (7)). Two nutrient types, oxygen and glucose, are defined. The x and y boundary for the nutrients are set to periodic to simulate a larger system. At the lower z boundary a Neumann boundary condition is applied simulating the impermeable ground and at the upper boundary a Dirichlet condition is applied simulating the connection to the bulk liquid. The initial condition and the inlet concentration are defined. The yield is set to 0.8 and the maintenance and decay rates are set to zero. The corresponding model parameters are depicted in Table 1. The microbe growth (see Line 5.2) is performed based on the Monod kinetics. The growth rate for the aerobic growth of the HET cells is calculated with the concentration of glucose S g and oxygen S O 2 as for currentNeighbour = 0 : numberOfNeighbours do 3: f

Parameter Symbol Value Unit
Glucose concentration (initial, boundary) S g,0 2 g/L Oxygen concentration (initial, boundary) First, the interaction between the cells is defined. Here, in contrast to the used contact force F c in NUFEB, a Lennard-Jones potential, which is already implemented in LAMMPS, is chosen. The standard 12/6 Lennard-Jones potential with a cut-off radius is computed according to where σ is the zero-crossing distance for the potential and the depth of the potential well. For the simulation exemplary values have been chosen. Next, for the interaction between the cells the experimental adhesion forces are imported from a table and considered for the withdrawal of attached cells. For the approach of the bacteria towards the wall a Lennard-Jones potential according to Equation (10) is defined. For the withdrawal of attached cells from the wall a table with the measured adhesion forces is imported. For the wall interaction the lower z plane is defined. At last, the planktonic movement is performed, where the z plane is specified for wall interaction and the maximum distance moved in one time step and the distance from the wall are set. A viscous force is implemented for stabilisation which is proportional to the velocity F i = −γ v i where γ is a damping coefficient.

Simulation
To obtain a first impression a sensitivity analysis is performed for which all parameters from Table 1 are independently added and subtracted by 10 % of the initial value. For the simulation one single HET cell is placed in the center of a domain and the biomass after 15.83 h is compared with the final biomass of the initial configuration. Next, the biofilm colonisation process is simulated. The measured force-distance curve from Figure 1 belonging to the species S. oralis is imported for the interaction of if type currentCell = planktonic and NumberOfNeighbours = 0 and position currentCell != wall then 3: perform random movement 4: else if type currentCell = planktonic and NumberofNeighbours != 0 and position currentCell != wall then 5: for currentNeighbour = 0 : numberOfNeighbours do 6: if type currentNeighbour = biofilm then 7: type currentCell = biofilm 8: end if 9: end for 10: if no neighbour = biofilm then 11: perform random movement 12: end if 13: else if type currentCell = planktonic and position currentCell = wall then 14: type currentCell = biofilm 15: end if 16 cells with the wall as the nutrient parameters have also been adapted for this species. Furthermore, the same table is also taken for the adhesion between the cells because no measured data between cells are available. The IbM parameters are taken from Table 1 except for the maximal growth rate which is doubled to 0.64/h. At the beginning of the simulation, 10000 planktonic cells are randomly placed in the domain. The initial diameter and density of the cells are set to 7.5 × 10 −7 m and 550 g/m 3 . Then, the attachment and growth is simulated for 16 h. The outcome of the cell distribution and the nutrient distribution are visualised in ParaView 5.9.0. (New York, USA). Furthermore, the development of the biomass and the cell numbers of the biofilm-associated and planktonic cells are compared. The simulation is performed three times to average out the randomness of the initial distribution of the planktonic cells as well as of the placement of mother and daughter cells.

Sensitivity analysis
In Figure 3 the results from the sensitivity analysis are depicted. It shows that the maximal growth rate and the initial/boundary concentrations as well as the halfvelocity constant of glucose have the largest influence on the final outcome of the biomass. The variations in positive and negative direction can be explained by referring to Equation (9). The diffusion coefficients do not have any influence on the outcome. This can be explained as the concentration of glucose and oxygen is large enough for the consumption of a single cell. The nutrients can diffuse faster than they are consumed by the cell. In a simulation with more cells the diffusion coefficient may have an influence on the outcome because the total nutrient consumption is larger.

Biofilm colonisation
The resulting colonisation and growth of the biofilm resembles the stages of the biofilm growth observed in nature. In Figure 4(a) the start of the simulation with 10000 randomly distributed planktonic cells can be seen. First, planktonic bacteria attach to the surface where the motile cells turn to immotile biofilm-associated cells (see Figure 4(b)). Here, they start to grow by consuming nutrients and small microcolonies are formed (see Figure 4(c)- Figure 4(d)). Next, these microcolonies continue to grow as cells multiply and the planktonic cell numbers decrease (see Figure 4(e)). Some of the microcolonies connect to form a larger biofilm (see Figure 4(f)) which also coincides with observations in nature [30]. The exponential increase in the numbers of living biofilm-associated cells and their biomass could be represented by the simulation (see Figure 5). The number of planktonic cells decreases as more and more biofilm-associated cells form. The average computation time of the simulation with one processor is 19 minutes with a standard deviation of 1 minute.
In summary, the model can be used to depict the qualitative adhesion of bacterial cells and the beginning of the biofilm growth. Several simulation trials (≥ 3) are needed to average out the randomness of the movement of the planktonic cells and the placement of the daughter cells after division.

Discussion
To improve the representation of the complex adhesion process, three new operations have been implemented   [3], so that it can be said they are plausible. Some commonly used forces in microbial simulations are available in NUFEB: a (wall) contact force, an EPS adhesion force and a drag force. The contact force is modelled via springs and dashpots and is only effective when cells overlap to eliminate the overlap during mechanical relaxation. The adhesion between the cells is modelled via a van der Waals force resulting only from the discrete EPS envelope or EPS particles and is neglected for cells that do not excrete EPS [15,24]. The parameters for these force calculations cannot be measured directly and have to be adjusted using calibration experiments. It is easier to rely directly on experimental data e.g. to use measured force-distance curves for the cell interaction which is pursued in this project. Instead of the contact force a Lennard-Jones potential is implemented for the cell-cell and cell-wall interaction which contains attractive and repulsive forces. Thus, the possibility to dissolve the overlaps of cells and the possibility of attraction between the bacteria and the wall is which are theories to explain the adhesion process considering the van der Waals and Coulomb forces as well as the hydrophobic/hydrophylic interactions, could be used [14]. The DLVO potential is already implemented in LAMMPS [27]. Nevertheless, both these potentials would have to be adapted by a calibration experiment to model the bacterial adhesion. It would be desirable to include measured force-distance curves for the approach of the cells as well but this cannot be measured yet. Wolff and Rudd [34] implemented tabulated potentials for the interaction between two particles so that the measured forces for the approach could easily be implemented in the future.
Furthermore, experimentally measured adhesion forces are included when a cell has previously attached to the wall and is retracted again and for the cell-cell interaction. This offers the possibility to rely directly on measured adhesion forces for the withdrawal of cells. A few simplifications have been assumed. The cell-wall and cell-cell interaction are not dependent on the size of the bacteria. In reality, one would assume that the adhesion forces depend on the size of the bacteria e.g. contact surface, radius or mass. Furthermore, the timedependency of the attachment process is not considered yet. In the IbM a cell has either attached or not and the attachment is simulated as instant where in reality it is time-dependent [2]. More and stronger bonds are formed over time with a limit reached after a distinct time. The reversibility of the adhesion is modelled as cells can be detached by large forces.
For the cell-cell interaction no distinction between approach and withdrawal is considered yet. The cells always consider the measured force-distance data when being withdrawn from other cells within a certain distance and not only if they have previously attached to another cell. For most cases the cells have attached to other cells before they are withdrawn e.g. a daughter cell which is placed next to its mother cell and therefore attached to it. Only for a few rare events two biofilmassociated cells have not attached before e.g. if two biofilm-associated cells are pushed towards each other by the fluid flow. The distance between the cells could be small enough for the adhesion forces to be considered without having attached to the other cell when they are separated again. The tabulated adhesion forces for the cell-cell interaction are only considered if no overlaps between cells occur. This ensures the establishment of a plausible starting position and in a next time step the adhesion forces can be considered. An external force could be large enough to separate the neighbouring cells without considering the adhesion force. This can be prevented by setting the maximal distance moved in one time step. Planktonic cells have been implemented by introducing a cell type which performs a step into a random direction with a random distance. They can be placed at a random location at the beginning of the simulation. Planktonic cells can become immotile if they are close to the wall or close to a biofilm-associated cell. Concerning the planktonic cells, it is assumed that every approach close enough to the wall leads to the attachment of the bacterium and they stay at this wall until an external force withdraws them. If the planktonic cells are in the vicinity of the biofilm, they turn into biofilm-associated cells. Biofilm-associated cells cannot turn to planktonic cells again. These cells remain biofilm-associated cells and are transported with the fluid flow. In nature, not all planktonic cells turn to biofilm-associated cells when they come close to the wall or other biofilm-associated cells. Furthermore, biofilm-associated cells can also return to their mobile lifestyle. Sweeney et al. [31] have included planktonic cells in the IbM solver iDynoMiCS. Therefore, they introduced a new agent which moves a random distance with a maximum distance at a random angle. If the bacterium is within a certain distance of a biofilm-associated cell, it switches its cell type. Furthermore, a chemotactic behaviour is included so that the planktonic cells can sense the surrounding nutrient gradient. They are repelled or attracted below or above a certain threshold to move in the designated direction.
The planktonic cells can also leave the biofilm which is regulated via a leaving probability that is influenced by local concentrations above the threshold. The IbM has been set up for the bacterial species of S. oralis which are primary colonisers in the oral cavity and provide attachment sites for secondary or late colonisers which can be pathogenic. Therefore, it is of interest to either prevent this initial attachment to the surface or the consecutive attachment of the pathogenic bacteria. All in all, the model provides qualitatively correct results for the attachment process and primal biofilm colonisation. It can be used to improve the understanding of the complex proceedings during the biofilm growth.

Future direction
Various extensions of the IbM to improve the representation of the bacterial adhesion and biofilm growth are reasonable.
Firstly, experimental data have to be generated for the validation of the model. This would include experimental data to adjust the parameters for the simulation e.g. the Lennard-Jones potential or the planktonic movement. Until now no external forces have been considered in the simulation. To finalise the validation of the implementation, the IbM should be extended to include external forces. Most of the time the adhesion of pathogenic bacteria or bacteria providing attachment sites for pathogenic bacteria should be prevented. For further analysis of this prevention different surface topologies could be implemented to discover how they influence the adhesion. Also different surface properties could be included by using force-distance curves obtained from the adhesion on different materials. Simulations could be used to analyse the influences on attachment for different surfaces and different bacterial species. In example, Palmer et. al [29] suggest that the adhesion might be favoured if the topography is of the same diameter as the cell diameter. An extension of the IbM to include a representation of EPS could be possible. So far in NUFEB EPS can be added as a discrete envelope around the cells which is secreted as an individual particle if a certain threshold is reached [24]. In contrast to discrete particles the EPS could also be modeled as continuous [16]. As the presence of EPS increases the bacterial adhesion, a dependence of the adhesion forces from the surrounding EPS concentration could be included [10]. Furthermore, the EPS secretion and beginning of the growth of the bacteria is mediated via autoinducers which could also be included in the simulation [32]. By introducing the production of autoinducers the EPS secretion and growth could begin above a certain threshold. In future simulations, the impact of fluid flow on the adhesion could be investigated. NUFEB offers the coupling with OpenFOAM so that a fluid simulation can be set up and a drag force, that acts on the bacterial cells, is calculated. In the oral cavity the influence of the salivary flow could be investigated and a constant fluid flow might be interesting in other areas e.g. for biofilm growth in industrial pipes leading to biocorrosion. Furthermore, a multispecies simulation could be performed, simulating e.g. the growth of a biofilm in the oral cavity. The used bacteria S. oralis are first colonisers in the oral cavity and thus an accessory pathogen. They provide attachment sites for the late pathogenic colonisers e.g. P. gingivalis. The corresponding force-distance curves between the cells have to be measured so that the attachment process of the multispecies biofilm can be analysed.
In the future, it may be possible to measure the forcedistance curves for the approach of bacteria to the wall or between bacteria. This would facilitate the usage of directly measured data and prevent the dependence on parameters that are calibrated with experiments for e.g. the Lennard-Jones potentials or the contact force. In addition, adhesion is assumed to be important throughout the biofilm development [16]. The model could be used to analyse the influence of adhesion throughout all stages of the biofilm development.

Conclusion
Bacterial biofilms can cause disadvantageous reactions which are linked to the adhesion of pathogenic bacteria on implants or the host tissue. Therefore, it is of interest to provide strategies to prevent the attachment of bacteria. The constructed Individual based Model is able to represent certain aspects of the complex bacterial adhesion. Measured adhesion forces can be imported for the interaction between the wall and between cells. Furthermore, the random movement of planktonic cells can be incorporated. It is demonstrated that the attachment of bacterial cells and the initial growth could be qualitatively modelled.
Code Availability: The commented source code and the input scripts are available under: 10.10.14464/gammas.v4i1.489.