A Thin-Film Lubrication Model for Biofilm Expansion Under Strong Adhesion

Understanding microbial biofilm growth is important to public health, because biofilms are a leading cause of persistent clinical infections. In this paper, we develop a thin-film model for microbial biofilm growth on a solid substratum to which it adheres strongly. We model biofilms as two-phase viscous fluid mixtures of living cells and extracellular fluid. The model tracks the movement, depletion, and uptake of nutrients explicitly, and incorporates cell proliferation via a nutrient-dependent source term. Notably, our thin-film reduction is two-dimensional and includes the vertical dependence of cell volume fraction. Numerical solutions show that this vertical dependence is weak for biologically-feasible parameters, reinforcing results from previous models in which this dependence was neglected. We exploit this weak dependence by writing and solving a simplified one-dimensional model that is computationally more efficient than the full model. We use both the one and two-dimensional models to predict how model parameters affect expansion speed and biofilm thickness. This analysis reveals that expansion speed depends on cell proliferation, nutrient availability, cell–cell adhesion on the upper surface, and slip on the biofilm–substratum interface. Our numerical solutions provide a means to qualitatively distinguish between the extensional flow and lubrication regimes, and quantitative predictions that can be tested in future experiments.


Introduction
Biofilms are sticky communities of micro-organisms that exist on surfaces. An estimated 2 80% of all microbes on Earth exist in biofilm colonies [1], and they have extensive effects on Given these assumptions, the general mass and momentum balance equations for our model are (2.1f) These equations differ from Tam et al. [12] only in that we allow nutrients in the biofilm to advect with both fluid phases, as opposed to assuming that the nutrients advect with the extracellular fluid phase only. In the mass balance equations, we adopt linear growth kinetics, which are the simplest forms that model cell proliferation proportional to cell volume fraction and nutrient concentration, and cell death proportional to local volume fraction. We also assume that dead cells immediately become part of the extracellular phase. The constants ψ n , ψ m , and ψ d in the source terms are the cell production rate, ECM production rate, and cell death rate, respectively. In the nutrient mass balance equations, the parameters D s and D b are the diffusivities of nutrients in the substratum and biofilm, respectively, and η is the nutrient consumption rate. We assume that nutrients disperse by diffusion in the substratum, and by both diffusion and advection inside the biofilm, where they can also be consumed by cells. To obtain the momentum balance equation (2.1f), we adopt a Stokes flow description and neglect inertial and body forces.
In axisymmetric cylindrical geometry, the relevant components of the stress tensor, σ, are (2.2) Due to cell proliferation and death, the stress components (2.2) retain divergence terms that 134 would otherwise vanish due to incompressibility. Invoking Stokes' hypothesis, we adopt 135 the standard coefficient −2µ/3 for these divergence terms [39,46]. We also neglect growth 136 pressure due to cell-cell contact, which was previously considered in similar models [51, 137 52]. Instead, we assume that microbes cannot respond actively to chemical or mechanical 138 cues from the environment. This is appropriate for yeasts because they are non-motile, 139 and also bacteria because they often lose swimming motility in biofilm environments [39]. 140 Instead of this growth pressure, we suggest that material incompressibility is sufficient to 141 drive expansion when cells proliferate.

142
To obtain general boundary conditions, we first assume that nutrients cannot pass through the base of the substratum or the biofilm-air interface. Nutrients can cross the biofilm-substratum interface, with a flux proportional to the local concentration difference across the interface. For the fluids, we assume that both phases exist in the biofilm only and cannot pass through the biofilm-substratum interface, on which we also impose a general slip condition. Finally, we impose the usual kinematic and zero tangential stress conditions on the free surface, and assume that free surface normal stress is proportional to mean local curvature. The complete boundary conditions are then t · (φ α σ ·n) = λφ α u ·t, u z = 0 on z = 0, (2.3c) where G is the initial nutrient concentration in the substratum. We then exploit the thin biofilm aspect ratio to obtain a simplified leading-order system of governing equations. Expanding variables in powers of ε 2 and applying the thin-film reduction, we obtain the leading-order governing equations (where variables are now expressed as dimensionless leading-order quantities) where the dimensionless parameters are (

Precursor Film Regularisation
Imposing a no-slip condition on the biofilm-substratum interface prevents biofilm expansion in the absence of suitable regularisation. This apparent paradox is commonly encountered in models involving a lubrication equation [41,46,53,54]. One method for dealing with moving contact lines is to introduce a precursor film [54]. This is an artificial thin layer of fluid with thickness b 1 existing ahead of the biofilm front. A physical interpretation of this is to represent the characteristic scale of surface roughness in the agar [53]. Following Ward and King [46], we adopt this precursor film to regularise the model. We assume that the precursor layer consists entirely of passive fluid (no cells), and that nutrient uptake or consumption does not occur in the precursor film. Mathematically, this involves modifying the model equations to extinguish the relevant terms wherever h ≤ h * , for some h * ≥ b [54]. The regularised model is then where θ denotes the Heaviside step function. The formulation (2.7) assumes that cells are 168 not present without sufficient biofilm thickness to support them, and that the biofilm cannot

172
Computing numerical solutions of the no-slip lubrication model then involves solving the 173 regularised system (2.7), subject to appropriate initial and boundary conditions.

Initial and Boundary Conditions
We obtain initial conditions for the regularised thin-film equations (2.7) by assuming that the substratum is initially filled uniformly with nutrients, and that no nutrients are yet present in the biofilm. For h(r, 0), we require a function such that there is a defined region with h = b ahead of the biofilm, and that relevant higher derivatives of h with respect to r are continuous throughout the entire domain. For this purpose, we modify the parabolic initial condition of Tam et al. [12] in the same way as Ward and King [46]. Finally, for the cell volume fraction φ n (r, z, 0), we choose a polynomial form such that φ n = 0 in the precursor film, and that the first derivatives of φ n with respect to r and z are continuous throughout the domain. Although this is not necessarily physical, it provides an initial volume fraction profile that varies with r and z. The initial conditions are then where H 0 is the initial height of the biofilm at its centre. The conditions (2.8) set the 176 initial contact line position to S(0) = 1, and satisfy h = b for r ≥ 1 and φ n = 0 for r ≥ 1.

177
For the boundary conditions, we assume radial symmetry at r = 0 in all variables.

178
Therefore, we obtain the conditions (2.9) 180 We also assume that the centre of the biofilm is fixed, and that fluid cannot cross

207
We solve the regularised two-dimensional lubrication model (2.7) numerically using the 208 Crank-Nicolson method. Rather than solve for φ n (r, z, t) directly, we instead introduce 209 the auxiliary variable which is the cumulative cell density through a vertical slice of the biofilm, measured from 212 the biofilm-substratum interface. We then haveφ n (r, t) = Φ n (r, h(r, t), t), and can recover 213 the solution for φ n using The numerical solution with initial conditions (2.8), boundary conditions (2.9)-(2.11), 225 and parameters listed in Table 2  radially, and attains a thicker shape to that seen in Ward and King [46]. As nutrients 234 deplete and cell proliferation decreases, surface tension forces facilitate slower radial spread,  We now exploit the observation in Figure 2.3 that cell volume fraction exhibits weak 247 dependence on z, and assume that φ n = φ n (r, t). This enables us to reduce the two-248 dimensional axisymmetric model (2.7) to a one-dimensional form, and eliminate u z .

249
Neglecting the two-dimensional structure has the major advantage of saving computational 250 time. After assuming that φ n is independent of z, equation (   The one-dimensional numerical solution using the parameters listed in Table 2

283
To investigate the effect of parameters on biofilm expansion speed and shape, we perform 284 a local sensitivity analysis. In each solution set, we vary one parameter from those in 285 Table 2.1, compute a numerical solution to t = 25, and calculate the biofilm radius and 286 thickness. To investigate the effect of cell proliferation rate, we vary the dimensional 287 parameter ψ n , which is otherwise scaled out of the model. We then update D, Pe, Ψ m , 288 and γ * in each simulation based on the value of ψ n , and compute solutions until T = ψ n /2. 289 Results for ψ n , Ψ m , and Ψ d then enable us to directly compare the effects on growth of 290 cell production rate, ECM production rate, and cell death rate respectively. Numerical  In the extensional flow model, we found that faster biomass production rates and increased 300 access to nutrients at the leading edge increased biofilm expansion speed [12]. Expansion 301 speed also depended on the initial condition, such that decreasing the initial biofilm height,     Local sensitivity analyses for biofilm radius at t = 25, with respect to net biomass production rates, initial biofilm height, surface tension coefficient, and slip parameter. We used the initial conditions (2.8) for all two-dimensional solutions, and replaced (2.8b) with (3.2) for one-dimensional solutions. Unless specified, all parameters used were those given in Table 2.1. Crosses represent results for the full two-dimensional model, and circles represent results for the simplified one-dimensional model.   Local sensitivity analyses for biofilm radius at t = 25, with respect to parameters that govern the movement, consumption, and uptake of nutrients. We used the initial conditions (2.8) for all two-dimensional solutions, and replaced (2.8b) with (3.2) for one-dimensional solutions. Unless specified, all parameters used were those given in Table 2 to the cells, which slows expansion. Conversely, increasing nutrient uptake rate, Q b , aids 367 cell production, because more nutrients become available for consumption.

368
In summary, expansion speed in strongly-adhesive biofilms depends on a combination 369 of cell proliferation (mediated by nutrient availability), surface tension forces, and slip.

370
Parameter changes that increase cell proliferation cause faster expansion, because they 371 increase the quantity of biomass. In contrast, increased surface tension forces, which 372 represents increased cell-cell adhesion on the upper surface of the biofilm, redistribute 373 biomass radially. This generates faster expansion, but also leads to thinner biofilms.
Similarly, the addition of slip on the biofilm-substratum interface promotes radial expansion 375 at the expense of vertical growth, giving rise to faster expansion. We reinforce these 376 conclusions in §4.2, where we investigate the effect of parameters on biofilm thickness.  Local sensitivity analyses for biofilm aspect ratio at t = 25, with respect to net biomass production rates, initial biofilm height, surface tension coefficient, and slip parameter. We used the initial conditions (2.8) for all two-dimensional solutions, and replaced (2.8b) with (3.2) for one-dimensional solutions. Unless specified, all parameters used were those given in Table 2.1. Crosses represent results for the full two-dimensional model, and circles represent results for the simplified one-dimensional model.  .4: Local sensitivity analyses for biofilm aspect ratio at t = 25, with respect to parameters that govern the movement, consumption, and uptake of nutrients. We used the initial conditions (2.8) for all two-dimensional solutions, and replaced (2.8b) with (3.2) for one-dimensional solutions. Unless specified, all parameters used were those given in Table 2.1. Crosses represent results for the full two-dimensional model, and circles represent results for the simplified one-dimensional model. Ψ m , Ψ d , and all parameters in Figure 4.4, do so because they affect the quantity of biomass. Increasing the quantity of biomass promotes both radial growth and thickening 388 of the biofilm, and therefore changes that increase expansion speed also increase biofilm 389 thickness. In contrast, the surface tension coefficient γ * and the slip parameter λ * affect the 390 distribution of biomass in the biofilm, and not its quantity. Since increasing γ * increased 391 the radial size, this occurs in conjunction with a reduction in thickness, as Figure 4.3e 392 shows. Similarly, decreasing the resistance to slip promotes radial expansion as opposed 393 to thickening, and therefore thickness increases as 1/λ * increases. advantage of the one-dimensional model is that it is computationally less expensive, and 415 potentially more amenable to analysis, than the two-dimensional model.