Coupled CFD‐DEM modeling to predict how EPS affects bacterial biofilm deformation, recovery and detachment under flow conditions

Abstract The deformation and detachment of bacterial biofilm are related to the structural and mechanical properties of the biofilm itself. Extracellular polymeric substances (EPS) play an important role on keeping the mechanical stability of biofilms. The understanding of biofilm mechanics and detachment can help to reveal biofilm survival mechanisms under fluid shear and provide insight about what flows might be needed to remove biofilm in a cleaning cycle or for a ship to remove biofilms. However, how the EPS may affect biofilm mechanics and its deformation in flow conditions remains elusive. To address this, a coupled computational fluid dynamic– discrete element method (CFD‐DEM) model was developed. The mechanisms of biofilm detachment, such as erosion and sloughing have been revealed by imposing hydrodynamic fluid flow at different velocities and loading rates. The model, which also allows adjustment of the proportion of different functional groups of microorganisms in the biofilm, enables the study of the contribution of EPS toward biofilm resistance to fluid shear stress. Furthermore, the stress–strain curves during biofilm deformation have been captured by loading and unloading fluid shear stress to study the viscoelastic properties of the biofilm. Our predicted emergent viscoelastic properties of biofilms were consistent with relevant experimental measurements.


| INTRODUCTION
Bacterial biofilms are initiated by reversible attachment of planktonic bacteria to a surface. Bacteria are then irreversibly attached to the surface and develop cell-cell cohesion. Matured biofilms are embedded in extracellular polymeric substances (EPS) which are produced by bacteria themselves (Flemming & Wingender, 2010).
The formation of biofilm helps bacteria to survive in harsh environments such as fluid flows (Banerjee et al., 2020). It was found that bacteria in biofilms are much more resistant to antibiotics than in planktonic state (Davies, 2003). Biofilms have dramatic impacts for a wide range of industries. For example, biofilms play an important role in bioremediation since they are able to convert toxic pollutants to harmless products (Singh et al., 2006;Yadav & Sanyal, 2019). Biofilms are also essential in wastewater treatment (Capdeville & Rols, 1992;Sehar & Naz, 2016). However, the accumulation of biofilms in industrial pipelines and drinking water systems may lead to biocorrosion (Abe et al., 2012;Klapper et al., 2002). Additionally, biofilms adhered to marine surfaces is an important trigger of accelerated biofouling (Antunes et al., 2019). The biofilms attached to the ship hull increase the frictional drag resulting in higher fuel consumption (de Carvalho, 2018). The emergence of biofilms allows pathogenic bacteria to survive in diverse environments (Tasneem et al., 2018). Besides, pathogen transmission is of concern to public health and can cause infection when the cells detach from the biofilm (Brindle et al., 2011). Therefore, a greater understanding of biofilm detachment in different hydrodynamic conditions may help to control the biofilm-related infection .
In biofilms, the EPS is a self-produced matrix which majorly consists of polysaccharides, extracellular DNA (eDNA) and protein (Erskine et al., 2018;Gloag et al., 2013;Yadav & Sanyal, 2019). It provides many functions, such as adhesion to surfaces and cohesion to maintain the mechanical stability of the biofilm system (Flemming et al., 2007). The production of EPS is essential during biofilm development since bacteria cells could be immobilized by EPS (Flemming & Wingender, 2010). EPS production can be responsively regulated, for example, it was found that EPS production could be affected by EPS biosynthetic genes (Ali et al., 2000;Song et al., 2018). Besides, mutant strains could cause the overproduction of EPS to help the biofilm position in the beneficial environment (Hibbing et al., 2010). All these could affect the EPS amount in biofilms. Biofilms may also increase the strength of the matrix by increasing EPS production when subjected to mechanical stresses at intermediate time scales (e.g., 1 h) (Shaw et al., 2004). Different biofilms can have different EPS and different mechanical properties (Houari et al., 2008;Klapper et al., 2002;Rupp et al., 2005;Stoodley et al., 1999;Vinogradov et al., 2004;Wloka et al., 2004). However, it is difficult to quantify EPS by microscopy or chemical analysis due to the complexity of the chemistry, as well as bias in extraction and purification techniques. Although EPS is complex, computational modeling can be simplified to represent its overall physical function rather than identify the individual polymer components, hence, gaining better understanding of the contribution of EPS production to biofilm mechanical properties.
In this study, a three-dimensional individual-based model (IbM) of biofilm was developed by coupling the computational fluid dynamics approach (CFD) with the discrete element method (DEM).
This model was implemented on NUFEB (https://github.com/nufeb) which is an open-source tool for individual-based modeling of microbial communities (Li et al., 2019). NUFEB integrates CFD-DEM solver SediFoam (https://github.com/xiaoh/sediFoam) which provides a flexible interface between large-scale atomic/molecular massively parallel simulator (LAMMPS) (Plimpton, 1995) and opensource field operation and manipulation (OpenFOAM) (Greenshields, 2017). The framework enabled us to describe the fluid induced biofilm deformation and detachment subjected to different flow velocities. In this study, we modeled a bacterial mutant that can produce the same type of EPS at different levels. Different EPS amounts were obtained by varying the relevant kinetic parameters in the model. We predicted the effect of EPS amount on the mechanical properties of biofilms and biofilm detachment.

| METHODOLOGY
The processes of biofilm growth and biofilm deformation were decoupled in this study, that is, fluid flow was applied to a pregrown biofilm. The pregrown biofilm was "grown" under static conditions for 5.3 days using the NUFEB individual based model which was described in Jayathilake, et al. (2017). The kinetic parameters for biofilm growth are provided in supporting information (Table S1). Only bacteria growth, division, and EPS production were considered in this study. Then the two-way coupling between the solid biofilm and computational fluid dynamic was adopted to investigate the deformation and detachment of biofilm under different hydrodynamic conditions. The simulation domain is displayed in Figure 1, with the pregrown biofilm positioned on the inlet side of the channel. The diameter of the involved particles is in the micrometry range (0.7-1.4 μm) based on the stochasticity of the biological system (Jayathilake, et al., 2017)

| Fluid-induced biofilm deformation and detachment
Experimental work has shown that EPS production in biofilms varied with bacterial strains and growth conditions (Costa et al., 2018;Danese et al., 2000). In this study, we achieved different EPS volume ratio (i.e., EPS volume divided by the volume of the biofilm) by changing the EPS growth yield coefficient in the modeling. The EPS growth yield coefficient was varied from 0.12 to 0.22 (g COD EPS /g COD S ) which corresponds to EPS/biofilm ratio of 20%-51% here. To investigate the biofilm deformation and detachment events, the biofilm with 46% EPS was subjected to inlet flow velocity between 0.1 and 0.4 m/s (Reynold number from 3.75 to 15, maximum wall shear stress from 10.7 to 42.7 Pa, shear rate from 2000 to 8000 s −1 ) for a duration of 40 ms. In the next simulation, the inlet fluid velocity was kept at 0.3 m/s to study the effect of EPS production on biofilm deformation and detachment. The detachment rate coefficient, which is defined as the ratio of the volume of detached biofilm clusters to the total volume of preformed biofilm, was calculated during the initial 14 ms (before biofilm washed away from the surface wall). Cluster detachment from the biofilm was defined as erosion if the particle number of the cluster was less than 1000 and sloughing if the particle number of the detached cluster exceeded 1000. The EPS amount, the mean and maximum heights, the roughness and porosity of different biofilms are summarized in Supporting Information (Table S2).

| Biofilm deformation-recovery test
The responses of biofilm to a rapid fluctuating shear stress were analysed immediately before biofilm failure. To save computational time, the fluid shear stress was applied to the biofilm for 3 ms (loading cycle) and then stopped immediately. Afterward, the biofilm was allowed to relax for 17 ms (unloading cycle). During the loading period, the fluid shear stress was increased by increasing the fluid velocity from 0 m/s at a constant acceleration. For the biofilm with 46% EPS, deformation-recovery tests were carried out by exposing the biofilm to the ramping flow with different accelerations: 20, 30, and 40 m/s 2 , respectively. Then the biofilms with 40% and 51% EPS were subject to the increasing fluid velocity at the acceleration of 40 m/s 2 to investigate the effect of EPS amount on mechanical response of biofilm. The shear strain in this simple shear test was defined as the angle change between the front edge of biofilm with the left channel wall ( Figure S1). The shear modulus was calculated as follows: where α is the shear strain, σ xz is the fluid induced shear stress on the biofilm which was computed globally by LAMMPS (Thompson et al., 2009). In this section, three planes (y = 5 μm, y = 15 μm, and y = 25 μm) were selected to measure the deformation angle thus to obtain the averaged shear strain ( Figure S2).

| Motion of bacterial and EPS agents
During biofilm deformation and detachment, the motion of bacterial cells and EPS agents are tracked by DEM on a Lagrangian framework: where v i ⃗ is the velocity of the particle i; m i is the particle mass; f c i , ⃗ is the contact force among collided particles (Xia et al., 2021), ⃗ is the fluid-particles interaction force.

| Locally averaged Navier-Stokes equations for fluids
The fluid flow is solved by locally averaged incompressible Navier-Stokes equation in which the fluid density ρ f is constant: ϵ s is solid volume fraction while ϵ f is fluid volume fraction which equals to (1 − ϵ s ). U s ⃗ and U f ⃗ are particle velocity and fluid velocity,

| Fluid-particle interaction
In this model, the fluid-particle interaction force f fp i , ⃗ consists of a drag force and lift force. For the particle i, the drag force model is expressed as (Sun et al., 2018): where V p i , is the volume of the particle i, U f i , ⃗ , and u p i , ⃗ are the fluid velocity and particle velocity, respectively. ϵ f i , is the fluid volume fraction while ϵ s i , is the solid volume fraction, β i is the drag correlation coefficient which is used to convert terminal velocity correlation to drag correlation (Syamlal et al., 1993).
In addition, the lift force on the particle i is calculated by the following formula (Sun et al., 2018;Van Rijn, 1984;Zhu et al., 2007): ⃗ is the curl of the flow velocity interpolated to the center of particle i.

| Cohesive force among particles
The cohesive force among the particles was computed by using the equation below (Israelachvili, 2011;Sun et al., 2018): where A is the cohesive strength, and s is the separation distance between the particle surface. A minimum separation distance s min was implemented when the separation distance between the two particles equals zero (s = 0). In this study, five different values of cohesive strength were used for the interactions of bacterial cells with bacterial cells, bacteria cells with EPS agents, bacteria cells with particle-wall, EPS agents with the particle-wall, EPS agents with EPS agents. Since EPS plays a significant role on binding the bacterial cells, the cohesive strength driven by EPS was assumed to be three orders of magnitude larger than that for bacteria (Fang et al., 2000).
The mechanical parameters of the simulations are listed in Table 1.

| EPS effect on biofilm deformation and detachment
To study how the EPS amount affected the deformation and detachment of biofilms, we examined a single fluid flow condition, 0.3 m/s sustained for a duration of 40 ms. For biofilms with a low amount of EPS (20% EPS), the biofilm clusters could easily detach from the biofilm matrix (erosion-dominated) at high frequency, accompanied by the escape of single bacterial cells ( Figure S6). This may be due to the limited EPS availability to immobilize the cells in biofilm (Flemming & Wingender, 2010).
The detachment frequency of biofilm decreased with the increase in EPS amount. Figure 4 shows the biofilms after being The local detachment, such as erosion and sloughing, could be significantly captured within the time of 14 ms. However, biofilm removal occurred after this period when the inlet fluid velocity was greater than 0.2 m/s. Therefore, the detachment rate coefficient, which is defined as the ratio of the volume of detached biofilm cluster to the total volume of biofilm per millisecond, was calculated during the initial detachment period (14 ms) and adopted to describe the biofilms detachment behavior under the range of hydrodynamic conditions. Each simulation was run for three replicates and the average results were calculated. As displayed in Figure 6a suggest that the resistance of the biofilm to the external fluid is largely attributable to the EPS amount. EPS is responsible for the mechanical stability of the biofilm due to its cohesive properties, therefore, the biofilms with a greater density of EPS components are predicted to be more stable when exposed to the fluid flow.

| Biofilm viscoelastic response during deformation-recovery test
Deformation of the biofilm with 46% EPS was monitored for 3 ms as the fluid velocity was incrementally increased from 0 m/s at a constant acceleration. Then the fluid flow was stopped, the biofilm was allowed to relax for 17 ms. The stress-strain curve was obtained from the loading and unloading cycle. Figure 7 shows the deforma-  Figure 8a shows the fluid-induced shear stress on biofilms overtime. It is evident that the higher flow acceleration resulted in higher shear stress imposed on biofilms (Figure 8a). This can lead to larger deformation (or shear strain) and deformation rate of biofilms (or strain rate), as seen in Figure 8b.
After the flow was removed at 3 ms, the fluid induced stress in the biofilm decreased rapidly (3-4 ms, Figure 8a), and some of the deformation (8%-10%) was immediately recovered attributable to a time-independent elastic response. Afterward (4-20 ms, Figure 8a), the stress decay slowed down dramatically and almost reached a plateau at the end of the recovery. A residual deformation (or strain) during biofilm relaxation was captured in each deformationrecovery test and increased with the maximum fluid velocity. Such a strain rate-dependent recovery was due to the nature of viscoelastic models adopted within the biofilms and is a common characteristic for viscoelastic materials (Capurro & Barberis, 2014;Chen et al., 2011). (or strain rate) due to the nature of viscoelastic effect. It is expected that the lower deformation rate leads to a smaller apparent shear modulus, which is also seen in this study. If we take into account different deformation rates of biofilms at different fluid velocities, it yields consistent equilibrium shear modulus and viscosity which are in the range of 3.8-4.8 Pa and 6.9-8.7 mPa*s, when using the Prony series viscoelastic model for curve fitting at different deformation rates (Chen & Lu, 2012).
To study the EPS effect on biofilm mechanics, we also focused on biofilms containing higher EPS amounts (40%, 46%, and 51%) subjected to the ramping fluid velocity at a constant acceleration of 40 m/s 2 . Biofilms with lowers EPS amount were not considered here since they could easily detach, thus the stress-strain curve would not be captured during deformation. The deformation of the biofilm with 40% EPS was greater compared to biofilm with 46% EPS, and both experienced similar shear stresses. This suggests that higher EPS resulted in the better resistance to the external fluid shear force. However, the shear strain of the 51% EPS biofilm was lowest although the maximum stress was almost 16% higher than for its counterpart biofilms. Taken together, the results suggest that biofilms with higher EPS might be stiffer, which agrees with what was found in Gloag et al. (2020). As seen in Figure 9c, the stress-strain curve was almost linear at very small strains (<0.1), which was also found in our experimental measurements for flow induced biofilm deformation of Bacillus subtilis ( Figure S8). When ignoring the deformation rate (or strain rate) effect, the apparent shear modulus at given loading conditions was 12.82 ± 2.03, 15.21 ± 1.94, and 17.18 ± 3.3 Pa for biofilms with 40% EPS, 46% EPS, and 51% EPS, respectively. When the deformation rates were accounted for, it yields the equilibrium shear modulus of 3.5 ± 0.7 Pa, 5.7 ± 1.7 Pa, and 5.6 ± 0.4 Pa for those three biofilms when using the Prony series viscoelastic model for curve fitting at different deformation rates (Chen & Lu, 2012). These are consistent with several biofilms such as Pseudomonas aeruginosa (Stoodley et al., 1999) and Staphylococcus aureus (Rupp et al., 2005). The reported corresponding viscosity for those three biofilms is 8.4 ± 4.5, 6.6 ± 1.6, 8.3 ± 3.2 mPa*s, respectively.
There is no evidence for correlation between the equilibrium shear modulus and viscosity of biofilms, as found in other recent studies (Houari et al., 2008;Klapper et al., 2002;Rupp et al., 2005;Safari et al., 2015;Stoodley et al., 1999;Vinogradov et al., 2004;Wloka et al., 2004). The predicted viscosity was the resultant of the mechanical interactions for bacteria-bacteria, bacteria-EPS and EPS-EPS.
In general, all the simulated biofilms exhibited some strain stiffening effect followed by strain softening at larger strains, which is The fluid induced (a) stress and (b) strain on biofilms changed with time and the corresponding averaged (c) stress-strain curve when the flow was applied and terminated (standard deviation did not give here for the high resolution). The biofilms with 40%, 46%, and 51% extracellular polymeric substance were selected. The fluid velocity was applied at an acceleration of 40 m/s 2 . The error bars represent standard deviations based on three replicates.
due to the viscoelastic properties and the change in the microstructure of biofilms during the deformation (Figure 9c). This is consistent with experimental measurements of biofilms at a wide range of strains (Jana et al., 2020). ANOVA test was used for statistical analysis (α = 0.05). The result (p < 0.05) suggests that the change of shear modulus with EPS amount has a statistically significant difference.
After the fluid was stopped, the stress on the biofilms decayed exponentially and the deformation recovered slowly which is a common feature for viscoelastic materials (Chen & Lu, 2012). The overall biofilm deformation recovery was 16%, 20%, and 22% for biofilms with 40%, 46%, and 51% EPS within the simulation period, respectively. The results suggest that the abundant presence of EPS in the deformed biofilms makes a significant contribution to their mechanical recovery, as the bacterial cells are much more loosely associated with other bacterial cells than EPS agents. Similar results were also verified in the experimental work, which released that the EPS is required to induce bacterial rearrangement during stress relaxation (Peterson et al., 2014).

| CONCLUSIONS
A CFD-DEM coupled model developed here has enabled us to predict biofilm deformation and detachment under varied hydrodynamic conditions. When the biofilm was exposed to a steady fluid proportional EPS (less than 32%), the biofilm easily detached from the surface and dispersion of individual cells was observed. In these cases, the limited amount of EPS was incapable of protecting the biofilm bacteria from shear stress. Biofilm detachment frequency decreased with the increase of EPS amount. The biofilms were stiffer at higher loading rate, which is a typical characteristic for viscoelastic materials. Such viscoelastic features of biofilms also led to the hysteresis loop (energy dissipation), which was predicted by the stress-strain curves in our simulations and experimental measurements (Rupp et al., 2005). The predicted shape of stress-strain curve during the flow induced deformation is similar to our measurements at a comparable flow velocity. Furthermore, we found that higher EPS amount led to a higher apparent shear modulus of the biofilm at given flow velocity. The equilibrium shear modulus was also higher when the EPS ratio was relatively high. In general, the predicted equilibrium shear modulus of the biofilms was in the range of 3.5-5.7 Pa, which is consistent with experimental measurements of P. aeruginosa and Staphylococcus aureus reported in literature (Rupp et al., 2005;Stoodley et al., 1999). The nonlinear stress-strain characteristics of biofilms at large strains predicted by the simulations were comparable with some key experimental findings by rheometer measurements (Jana et al., 2020).

AUTHOR CONTRIBUTIONS
Yuqing Xia, Pahala G. Jayathilake, Paul Stoodley, and Jinju Chen designed the research. Yuqing Xia performed the simulation work and acquired the data. Yuqing Xia and Jinju Chen did data analysis. Pahala G. Jayathilake, and Bowen Li contributed to data analysis. Yuqing Xia, Pahala G. Jayathilake, and Jinju Chen prepared the original draft. Jinju Chen and Pahala G. Jayathilake provided the overall guidance of the work. All the authors contributed to the writing of the manuscript, revised it, and approved the final version.