Using particle tracking, particle positions in each frame for the acquired images were obtained. The image acquisition frequency yielded particle trajectories or displacement as a function of time. Individual particle creep curves, giving material strain versus time, were plotted by converting the displacement curves into creep compliance curves (total compliance as a function of time J(t) ), knowing the probe deflection d(t) and the applied force F , according to the compliance formula (Schnurr, Gittes, MacKintosh, & Schmidt, 1997): , which gives the relation between material strain and applied stress for a probe of radius R (the microbead) embedded in an incompressible viscoelastic medium. The compliance (J(t) ) of the material is the inverse of elastic modulus and equal to the ratio of strain to stress, i.e., the deformation for a given stress. Thus, a higher compliance indicates a more flexible material. The absolute force from the magnetic field acting on the magnetic particles was determined by a purely viscous mixture of 99.5% glycerol with a concentration of 2\(\times\)10⁸ particles/mL (see SI).

To study the impacts of mechanical heterogeneity on biofilm deformation, a two-dimensional (2D) continuum biofilm model was implemented using COMSOL Multiphysics (COMSOL v4.4, Comsol Inc, Burlington, MA) with finite element method. A fluid-structure interaction (FSI) model using the arbitrary Lagrangian-Eulerian (ALE) approach was applied to simulate the deformation of biofilms under fluid flow. The biofilm was considered either a purely linear elastic material with elastic compliance varying with depth in the biofilm, as per the experimental data, or a homogenous material, using the average of the experimental data. The system included a solid biofilm subdomain, and a bulk liquid subdomain. Two arbitrary biofilm structures were considered: a biofilm colony and a mushroom-like tower. The spatial distribution of mechanical properties was set in the biofilm structure using the experimental data obtained for the growth under 0.1 mL/h flow rate, DO saturation, and no Ca²⁺ addition, with an elastic modulus (E) range of 0.005-0.05 Pa, as calculated from the compliance (J ) determinations. The averaged elastic modulus of the complete heterogeneous biofilm structure was computed over the biofilm domain and then implemented into the homogeneous model. Time-dependent simulation was performed for 2-3 seconds until obtaining a steady-state solution. Biofilm deformations at 1.5 s were reported. A detailed description of the FSI model is provided in the SI.

We used magnetic microbead actuation together with confocal fluorescence microscopy to assess the spatial distribution of mechanical properties in GFP-tagged P. aeruginosa PAO1 biofilms. Creep curves showed the type of mechanical response ranging from elastic to viscoelastic. Non-magnetic beads did not show significant motion in the respective control at distances greater than 5 μm. In the local viscoelastic environments, the deformation observed for the loading time was dominated by elastic compliance, i.e., the viscous part of the deformation was lower that the elastic initial deformation. Also, about 20 to 30% of the particles (depending on the condition) showed purely elastic response. Thus, the analysis was focused on the spatial distribution of the elastic compliance J for each biofilm growth condition. From J , the elastic modulus (E) can be obtained as its reciprocal and applied for mathematical simulations.

The effect of flow-induced shear stress on the mechanical properties was studied by growing biofilms under different flow rates. The distribution of compliance values in different directions is shown in Figure 2, for biofilms grown under 0.1, 1 and 5 mL/h (Re = 0.28, 2.8 and 13.9), respectively. Biofilms grown at the higher flow rates of 1 and 5 mL/h were thinner (approximately 80–150 μm) than at 0.1 mL/h (approximately 150–250 μm). The spatial distribution of compliance revealed significant variability in mechanical properties. Higher heterogeneity was observed for the lower flow rates, with compliance values differing by as much as two orders of magnitude within the same biofilm (Figure 2). The distribution of compliance showed a significant correlation with biofilm depth (Z), some dependence across the capillary width (X), and was variable but not stratified across the capillary length (Y).

The lowest flow rate, 0.1 mL/h (Re=0.28), produced the greatest vertical heterogeneity of compliance, with more elastic areas at the biofilm exterior and more rigid areas at the base. For the 1 mL/h flow rate (Re=2.8), a more evenly layered structure was observed. The biofilm grown at 5 mL/h (Re=13.9) showed high rigidity throughout, with a mean compliance value as low as 0.02 m²/N (E=50 Pa). Increasing the flow shear stress decreased the average vertical compliance from values of 2 m²/N to 0.02 m²/N (elastic moduli range of 0.5–50 Pa), producing more rigid structures. Similar results were observed in E. colibiofilms in a previous study using the same experimental technique (Galy et al., 2012). These results were observed for low shear stresses during growth (Re from 0.28–13.9 corresponds to a shear stress of 0.22–11 mPa). Thus, the shear stress during growth was not expected to cause deformation, i.e., the elastic modulus is high compared to the wall shear stress in the channel. This suggests that the observed biofilm stiffening may be explained by a biological adaptation to low shear rather than a physical remodeling due to shear stress (Galy et al., 2012).

The compliances of biofilm basal layers (0-20 μm) were almost the same for the different flow rates, despite greater differences in the outer layers. This suggests the mechanical properties near the biofilm base were not significantly affected by shear. This may be because the outer layer of the biofilm is more affected by the bulk conditions (Even et al., 2017). It also may be due to the lower biological activity at the base, leading to lower net growth rates compared to the upper layer. The deeper layers of biofilms are known to consolidate and became more dense over time (Laspidou & Rittmann, 2004a).

The compliance across the capillary width (x-direction) showed a lower value near the center of the capillary channel. Result showing increasing stiffness in the center of biofilms have been reported by others (Karampatzakis et al., 2017). The variation of compliance along the capillary width may be caused by the shear force distribution. Stiffer biofilms in the middle of the flow channel correlate with the higher shear stress condition, while lower shear near the edges led to a softer biofilm. This is consistent with the findings from Thomen et al. (2017) who found that the development of biofilm had a strong relationship with the spatial distribution of shear stress.

The non-uniformity of biofilm mechanical properties can be correlated with the differences in concentrations of biofilm constituents, due to changes in EPS density or water content (Wilking et al., 2011). For instance, different EPS polysaccharides, such as Psl and Pel (Yang et al., 2011), produced by biofilm microorganisms such as P. aeruginosa can either increase the elasticity or increase the viscosity of biofilms, respectively (Chew et al., 2014; Friedman & Kolter, 2004). Thomen et al. (2017) suggested that surrounding flows may partially wash out extracellular molecules secreted by cells. Based on the elastic results in our study, it is possible that shear flow removed bacteria cells instead of extracellular molecules, leading a higher stiffness under high shear stress. Although relatively low shear was used in the experiments, it is likely that continuous and higher shear stresses would result in an increasing production of EPS components such as Psl, increasing the strength of the matrix. Previous work also found an influence of hydrodynamics on the levels of quorum sensing molecules (Timp et al., 2009). Quorum sensing impacts gene expression and often promotes and regulation biofilm formation and EPS production. Therefore, it is also possible that shear stress changed the biofilm elasticity due to its impacts on quorum sensing.

The effects DO on biofilm mechanical properties were studied for two flow rates, 0.1 mL/h (Re=0.28) and 5 mL/h (Re=13.9), using the methodology previously described. Compared to the biofilms obtained under DO saturation (8 mg/L), biofilms grown at a DO of 1 mg/L showed similar thickness but were more compact (less rough) (Figure 3). The vertical distribution of compliance values for these biofilms after 5 days of growth is shown in Figure 3.

At equal flow rates, biofilms displayed higher compliance values at lower DO levels. Similar to biofilms at high DO, the biofilms at low DO produced a wider range of compliance values flow for the lowest flow rates (Re=0.28). The outer biofilm layers had a higher compliance than the basal layers (Figure 3). Interestingly, the highest flow rate (Re=13.9) also produced a layered biofilm, not as stiff and uniform as the one with high DO. The average compliance at Re=13.9, around 0.6 m²/N, was lower than for Re=0.28, around 1.2 m²/N. These results show the same trend discussed above, that increased shear stress leads to increased biofilm stiffness.

Higher compliances were observed in the superficial biofilm under low DO conditions. Previous research revealed a correlation between DO and growth rate and EPS production (Applegate & Bryers, 1991; Laspidou & Rittmann, 2004b). Laspidou and Rittmann (2004) used modeling to infer that consolidation in the bottom layers resulted in higher density and lower porosity. Lower production of polysaccharides was also associated with oxygen-limited conditions (Ahimou et al., 2007; Applegate & Bryers, 1991). Under different bulk DO conditions, oxygen depletion with depth caused a variation of cell density and EPS production, which increased the heterogeneity of mechanical properties. As mentioned above, higher compliance in the superficial biofilm at low DO may be caused by low growth rate of cells. Furthermore, several studies have demonstrated that the oxygen gradients in biofilms are determined by the biological activity (Bridier et al., 2015; Xu et al., 1998). Thus, biofilm properties in the bottom layer are not strongly correlated with the oxygen concentrations in the bulk. The consistency of compliance values in the bottom layer was likely due to the lack of growing cells and anoxic conditions, regardless of the DO in the bulk liquid.

The compliance profile for biofilms grown at the lowest growth flow rate (Re=0.28) was determined for presence of Ca²⁺, as shown in Figure 4. In both cases the compliance values decreased, especially in the DO-saturated condition, where instead of the widely distributed compliance observed without Ca²⁺ addition, the results showed a reduced spread in values and a more homogeneous compliance profile. The change in heterogeneity can be correlated with the chemical effects of Ca²⁺, which form a bridge between negatively charged moieties in the EPS (Flemming & Wingender, 2001). Ca²⁺ has been shown to promote ordered protein helices, and has a strong affinity for metal ions (Sehar et al., 2016; Sutherland, 2001). The change of compliance was due to the Ca²⁺, which promoted a stronger and stiffer biofilm matrix. Previous studies showed that Ca²⁺ could stimulate the development of thick and compact biofilms, with increased mechanical stability and elastic modulus (Ahimou et al., 2007; Körstgens et al., 2001; Shen et al., 2018). The results suggest that the development of stiffer structures was likely to be a consequence of high degrees of cross-linking.

The presence of Ca²⁺ led to much thicker biofilms when the bulk DO was high, and slightly thicker biofilms when the bulk DO was low. Sehar et al. (2016) concluded that the addition of Ca²⁺ could increase the cell density and thicken the biofilm. Thus, the formation of mechanically more stable biofilms in the presence of Ca²⁺ may be due to the strong cationic bridging between bacterial cells and EPS polymers. This is consistent with results obtained with other multivalent cations that play the same role as Ca²⁺ , such as Cu²⁺, Mg²⁺, and Fe³⁺ (Beech & Sunner, 2004).

Past research has mainly assessed the mechanical properties of biofilms with macroscale tools (e.g., shear rheometry), characterizing an entire biofilm with one reading. Other studies have characterized only the outer layers of the biofilm (e.g., microindentation). Also, few modeling studies have considered mechanical heterogeneity. To illustrate the potential impacts of mechanical heterogeneity on biofilm behavior (i.e., biofilm deformation), as well as the potential differences in assuming biofilms have homogeneous elastic properties, we used a 2D FSI model to simulate the mechanical responses of two arbitrary biofilm morphologies. We arbitrarily assigned the experimental observed z-variations of mechanical properties from this study. Two biofilm morphologies were considered: a biofilm cluster or colony, and mushroom-like biofilm tower (Figure 5a and 5b).

The biofilm structures were modeled as purely linear-elastic solids, i.e., without viscous behavior. The initial morphologies were arbitrary, and it was assumed that the timescales were short enough that morphology changes due to growth or decay could be neglected. Thus, the model captures the short-term behavior when exposed to fluid flow. For each case, the deformations were compared for uniform mechanical properties and properties varying in the z direction per one of our experimental results.

Figure 5a and 5b shows 2D plots of the fluid velocities, von Mises stresses in the biofilm, and initial and final biofilm positions. The velocity of biofilm colony simulation (Re=250) was two orders of magnitude higher than that of biofilm tower simulation (Re=1) since the tower shape is more easily deformed. The original and deformed biofilm structures are shown as white outlines and colored surfaces in Figure 5a and 5b.

To evaluate the impact of mechanical heterogeneity on biofilm deformation, the displacement of the biofilm colony was plotted for as a function of biofilm depth (z direction) (Figure 5a and c). For the biofilm colony (Figure 5c), the heterogeneous biofilm suffered larger deformation. The more deformable top layer in the heterogeneous biofilm showed three times greater deformation (55 µm at the top) than the homogeneous biofilm colony (20 µm at the top). This is a 64% increase in deformation. For the mushroom-like biofilm tower (Figures 5b and d), results showed larger deformations (70 µm at top) for the homogeneous biofilm than for the heterogeneous biofilm (55 µm at top), a difference of 22.8%.

The modeling simulations show that, even when the spatial distribution of elastic parameters is the same, the effect of mechanical heterogeneity on overall deformations can vary depending on the biofilm morphology. The assumption of homogeneous mechanical properties can lead to significant differences in deformation predictions.

Possible mechanisms and implications of biofilm mechanical heterogeneity

Our experimental results showed that P. aeruginosa biofilms have a heterogeneous distribution of mechanical properties, with a wide range of values. Other researchers have explored the mechanical non-uniformity both experimentally (Birjiniuk et al., 2014; Cao et al., 2016; Galy et al., 2012; Karampatzakis et al., 2017) and theoretically (Bridier et al., 2015; Even et al., 2017; Laspidou et al., 2014). The spatial distribution of biofilm mechanical properties was also found to vary with environmental conditions. The overall impacts of chemical, biological, and physical factors from the environment and biofilm itself explain the mechanical heterogeneity of the biofilm.

Stewart and Franklin (2008) concluded that mechanisms of biofilm heterogeneity including, but were not limited to, microscale chemical gradients, adaptation to local environmental conditions, stochastic gene expression, and genotypic variation. Another possible mechanism is that the growth of biofilms tended to be in the direction of minimal mechanical resistance, which lead to consistent gradients of stiffness, oxygen concentrations, nutrient concentrations, and growth rates (Even et al., 2017).

Biofilm heterogeneity has important practical implications, especially for biofilm control. Past research on the mechanical properties of biofilms show a wide range of results, possibly due to the wide variety of testing techniques at different length scales. In addition, spatial heterogeneity of mechanical properties, and differences in biofilm morphologies and structures, play key roles on the deformation of biofilms under fluid flow (Even et al., 2017; Trejo et al., 2013). Our modeling results show that the spatial heterogeneity of biofilm mechanical properties could lead to significant differences of the deformation, even with the same averaged value. For this reason, previous studies (Picioreanu et al., 2018; Stoodley et al., 1999) which back-calculated homogeneous mechanical properties using the experimental deformation should be viewed with caution.

Our modeling results show that the variability of mechanical properties can have different effects, depending on the biofilm morphology, spatial distribution of mechanical properties, and hydrodynamic profiles. Surface erosion and biofilm sloughing may occur more easily and frequently on more non-uniform biofilm colonies and more uniform mushroom-like biofilms.

In water and wastewater engineering, biofilms can play both beneficial and detrimental roles. For example, the accumulation of biofilm in biotreatment processes is critical to their good performance, while development of biofilm fouling layers in membrane filtration systems can greatly increase energy requirements. Thus, promoting accumulation of beneficial biofilms and promoting detachment of detrimental biofilms is an important practical need. Yet there are few tools available to predict biofilm behavior. This research provides both results on the variability of mechanical properties for a P. aeruginiosabiofilm, and also a tool to assess its impact on biofilm deformation.

In order to affect the mechanical stability of biofilms, the fluid flow rates and shear, DO concentrations, and concentration of divalent cations are important variables. Our results also suggest that biofilm initially may be more flexible and easy to detach at the base, but then become stiffer at the base as the biofilm matures. Thus, lighter and more frequent biofilm removal treatments may be more effective than ones that are less frequent and more intense. For longer term assessments of biofilm behavior, models should include growth and decay processes, and also the viscoelastic behavior of biofilms.

Magnetic tweezers were used to characterize the spatial distribution of elastic properties of P. aeruginosa biofilms for a range of conditions. Spatial heterogeneity was observed in all three-dimensions, and environmental conditions had a significant impact on the spatial distribution. Higher shear resulted in a stiffer and more uniform distribution, possibly due to mechanical adaptation. Also, stiffer biofilms in the center of the flow cell were correlated with greater hydrodynamic shear. Lower bulk DO led to a more heterogeneous biofilm, probably due to the greater variation of biological activity with depth in the biofilm. The addition of Ca²⁺ in bulk liquid increased the average stiffness and resulted in more uniform biofilms. Further research should address the viscoelastic behavior, and the combination of elastic deformation, viscous deformation, and growth, as all of these may significantly impact biofilm deformation and morphology in the long term.

Using mathematical model for two different hypothetical biofilm morphologies under fluid flow, it was shown that the spatial heterogeneity of mechanical properties can lead to significant differences in biofilm deformation. This demonstrates that biofilm mechanical heterogeneity should be considered when predicting biofilm deformation. Conversely, it also should be taken into consideration when using biofilm deformation to infer its mechanical properties.