Nutrient Source Location and Biolm Viscoelasticity Control Biolm Growth Rate, Migration Rate, and Morphology in Shear Flow

We present a numerical model to simulate the growth and deformation of a viscoelastic bioﬁlm in shear ﬂow under different nutrient conditions. The mechanical interaction between the bioﬁlm and the ﬂuid is computed using the Immersed Boundary Method with viscoelastic parameters determined a priori from measurements reported in the literature. Bioﬁlm growth occurs at the bioﬁlm-ﬂuid interface by a stochastic rule that depends on the local nutrient concentration. We compare the growth, migration, and morphology of viscoelastic bioﬁlms with a common relaxation time of 18 minutes over the range of elastic moduli 10-1000 Pa in different nearby nutrient source conﬁgurations. Simulations with shear ﬂow and an upstream or a downstream nutrient source indicate that soft bioﬁlms grow more if nutrients are downstream and stiff bioﬁlms grow more if nutrients are upstream. Also, soft bioﬁlms migrate faster than stiff bioﬁlms toward a downstream nutrient source, and although stiff bioﬁlms migrate toward an upstream nutrient source, soft bioﬁlms do not. Simulations without nutrients show that on the time scale of several hours, soft bioﬁlms develop irregular structures at the bioﬁlm-ﬂuid interface, but stiff bioﬁlms deform little. Our results agree with the biophysical principle that bioﬁlms can adapt to their mechanical and chemical environment by modulating their viscoelastic properties. We also compare the behavior of a purely elastic bioﬁlm to a viscoelastic bioﬁlm with the same elastic modulus of 50 Pa. We ﬁnd the elastic bioﬁlm underestimates growth rates and downstream migration rates if the nutrient source is downstream, and it overestimates growth rates and upstream migration rates if the nutrient source is upstream. Future modeling can use our comparison to identify errors that can occur by simulating bioﬁlms as purely elastic structures.


Introduction
Biofilms are ubiquitous in nature but could be detrimental in industrial and medical applications. When planktonic microorganisms attach to a solid surface immersed in a fluid flow, they become sessile and develop into biofilms 1 . Corrosion in pipelines due to biofilm formation is costly to the oil and gas industry, water utilities, and power plants 2 . In water transported to consumers through pipelines, biofilms growing on pipe walls release bacterial cells that deteriorate water quality and alter the drinking water microbiome 3 . The metal surfaces of vascular stents used to increase the diameter of arteries and improve blood flow in medical patients 4 are susceptible to biofilm formation, which could lead to inflammation, necrosis, and vessel rupture 5 . In the dairy industry, biofilms that form on the stainless steel surfaces of milk storage tanks and process lines increase the risk of microbial contamination of processed dairy products 6 . In all these industrial settings, biofilms grow in a fluid environment under different nutrient supply and flow conditions.
Biofilms are composed of microorganisms bound together by a matrix of polymers, where the polymers are secreted by the cells 7,8 . The matrix responds to mechanical stress by exhibiting elastic strain and viscous creep [9][10][11][12] and its chemical composition facilitates nutrient capture 13 . By modulating their viscoelasticity and chemical composition through genetic expression, biofilms are capable of adapting to their mechanical and chemical environments [14][15][16][17] . For example, biofilms are able to resist mechanical clearance in the lungs of cystic fibrosis patients by adjusting their polysaccharide production 18 and resist the immune response in chronic wound patients by producing rhamnolipids from quorum sensing cells 19 . Numerical models that incorporate realistic biofilm properties have the potential to make quantitative predictions about the mechanical and chemical behavior of biofilms in complex environmental conditions. Such predictions could suggest ways to mitigate biofilm growth in industrial and medical settings.
Mathematical models with different levels of complexity have been developed to simulate biofilm growth and deformation in fluid environments (for comprehensive reviews, see Refs. [20][21][22] and references therein). The models can be generally classified into two groups 22 : those that simulate deformation, and those that simulate growth. The Immersed Boundary Method 23 (IBM) has been a popular method of simulating biofilm deformation because it easily captures mechanical interactions between immersed structures and fluids. Some of the IBM models represented biofilms as purely elastic structures of nodes connected by Stokes elements 24,25 . Other models accounted for biofilm viscosity by changing the fluid viscosity near elastic biofilms 26,27 or by representing biofilms as nodes connected by viscoelastic Kelvin-Voigt elements 28 . However, none of these IBM models included growth. Furthermore, those that included viscosity did not exhibit a constant strain rate at long times, hence, they could not reproduce the viscoelastic relaxation times measured in experiments 9,10 . Another set of models used discrete stochastic rules to simulate growing biofilms as they consume nutrients from their environment [29][30][31][32][33][34][35][36][37][38] . These models, however, did not include biofilm deformation. For large-scale simulations of biofilm processes, continuum models are utilized to simulate the fluid flow, nutrient distribution, and the spatial expansion of biomass. Examples include phase-field approaches that simulate the growth of non-elastic viscous biofilms 39 or the deformation of viscoelastic biofilms that do not grow 40 . Even though these continuum models can produce results that match elaborate experiments quantitatively, their formulations tend to be complicated and computationally challenging 22 . Hybrid models that combine discrete biofilm structures and stochastic growth schemes with continuum formulations of nutrient and fluid dynamics can be easy-to-implement and computationally efficient.
In this paper, we present a novel hybrid biofilm model to study how the mechanical properties of biofilms affect their interactions with different nutrient supply and fluid flow conditions. The biofilm is represented as a structure of nodes connected by viscoelastic Maxwell elements with viscoelastic parameter values determined a priori from experimental measurements of shear modulus and relaxation time 10,41 . The 2D simulation is based on the IBM and utilizes the IBM solver in the open source IB2d software [42][43][44] . An advection-diffusion-reaction equation, coupled with the IBM, governs the nutrient concentration in the fluid. As the biofilm consumes nutrients, new nodes and visocleastic elements are added to the biofilm-fluid interface by a stochastic rule that is based on local concentration values. Our model can be used to predict the growth rate, migration rate, and morphology of viscoelastic biofilms that result from the interplay between viscoelastic deformation from fluid shear stress and growth from nutrient consumption over several hours.
Our paper is organized into three sections. In the Methods section, we discuss our approach to IBM that enables an a priori determination of parameter values, and we describe our biofilm growth and reaction algorithm. In the Results section, we present results from simulations with biofilms in shear flow with different flow velocities and nutrient source configurations. First we compare the results between a purely elastic biofilm and a viscoelastic biofilm with the same elastic modulus. Then we present the results for viscoelastic biofilms with a common relaxation time of 18 mintues 10 over the range of elastic moduli 10-1000 Pa 10,41 . In the Discussion we summarize the significance of our methods for future modeling work and interpret our results in a biophysical context.

Numerical Model
Biofilms immersed in fluid and attached to a surface deform in response to fluid flow as they grow by consuming nutrients from the fluid environment. The deformation response of the biofilm can be described by bulk viscoelastic material properties, and the nutrient response by spatiotemporal chemical kinetics. Here, we present a biofilom model to simulate the growth of a viscoelastic biofilm attached to the wall of a microchannel with fluid flowing through the channel and a nutrient source located along the wall.
We use the Immersed Boundary Method (IBM) 45 to simulate the hydrodynamic interactions between water and a growing viscoelastic biofilm, with viscoelastic simulation parameters determined a priori from experimental data. The fluid is defined in fixed Eulerian Cartesian coordinates x = (x, y) and coupled to the viscoelastic biofilm dynamics defined in variable Lagrangian Cartesian coordinates X(t) = (X(t),Y (t)). The Lagrangian coordinates of the biofilm are treated as nodes in a triangular mesh where the connections represent viscoeslastic elements that exert forces F * (X,t) on the nodes. In IBM, the forces exerted on the Lagrangian nodes are transformed into body forces f(x,t) exerted on the surrounding fluid particles in the Eulerian coordinates. The fluid attains velocity u(x,t) as a result of its interactions with the immersed structure and other external forces. The fluid velocity is then interpolated to the Lagrangian nodes of the biofilm giving them velocities U(X,t) that match the local fluid velocities so that the no-slip boundary condition is imposed at the immersed structure-fluid interface. Even though the fluid mechanics in our simulations are at low Reynolds number, and are thus dominated by viscous effects, we prefer to use the IBM to solve the Navier-Stokes equation so that we can simulate biofilms at high Reynolds numbers in the future.
The dynamics of the concentration c(x,t) of a nutrient affecting biofilm growth is governed by an advection-diffusionreaction equation solved in the Eulerian coordinates with a reaction term r(x,t) describing the nutrient uptake by the biofilm. As the biofilm consumes nutrients from the surrounding fluid environment, it grows by adding nodes and triangles to the Lagrangian structure via a stochastic rule based on local growth rates that are determined by local nutrient concentrations. In the following subsections, we outline the details of our approach.

Immersed Boundary Method
The motion of an incompressible Newtonian fluid that interacts with immersed structures is governed by the Navier Stokes equations where u is the fluid velocity, p is the pressure, ρ is the fluid density, µ is the dynamic viscosity, x is position in the Eulerian coordinates, and t is time. The net body force exerted on the fluid has a contribution f from the interaction between the immersed structure and the fluid and a contribution f ext arising from all other external forces.
The density of forces acting on the biofilm due to internal stresses, described in Lagrangian coordinates as F(X,t) at a position X, is expressed in Eulerian coordinates through the equation where the integral is taken over the Lagrangian domain L . The integral in Eq. 3 with the Dirac delta function kernal δ (x − X) transforms the force density distributed on the immersed Lagrangian structure into body forces that are distributed on the surrounding fluid in the Eulerian domain. When the no-slip boundary condition is enforced at the immersed structure-fluid interface, the immersed structure must move at the local fluid velocity. In this case, the velocity of the immersed structure is given by where the integral is taken over the Eulerian fluid domain E . The integral transform in Eq. 4 interpolates the velocity of the fluid in the Eulerian domain to the immersed structure in the Lagrangian domain.
We use the open-source software IB2d 42-44 to solve Eqs. 1, 2, 3 and 4 with periodic boundary conditions imposed at the inlet and outlet of the fluid domain. The Eulerian domain is discretized into a set of points x i j = (x i , y j ) defined on a square grid with grid spacing h. The Lagrangian structure is discretized into a set of points X m (t) = (X m (t),Y m (t)) where the distance between points may change at each time step. Eqs. 3 and 4 are implemented numerically by using the regularized delta function 45 which has compact support within a square of side length 2h.

Viscoelastic Biofilm-Wall Structure
The biofilm in our model consists of set of Lagrangian nodes connected as a triangular mesh. To model the mechanics of the biofilm, the connections between nodes are treated as viscoelastic elements. The elements may be either a Stokes element, which is a linear spring that obeys Hooke's law, or a Maxwell element, which consists of a Hookean spring and a dashpot in series. The force acting on the Lagrangian node X m due to the element connecting it to Lagrangian node X m ′ is given by where E is the elastic (Young's) modulus of the biofilm and d mm ′ (t) is the resting length of the element. If the elements are purely elastic springs, the resting length, d mm ′ (t), in Eq. 7 is set equal to the initial resting length, d mm ′ (0).

3/20
For viscoelastic elements, the resting length of the element, d mm ′ (t), varies in time as the dashpot stretches or compresses. The rate of change of the resting length,ḋ mm (t), is proportional to the magnitude of the spring force 46 and is given bẏ where η is the dashpot viscosity. In the limit η → ∞ at finite E,ḋ mm ′ (t) → 0, thus the element connecting X m and X m ′ becomes a purely elastic spring. The walls of the simulation channel are each composed of a linear chain of Lagrangian nodes of initially equal spacing. Neighboring nodes are connected by stiff springs and each node is tethered to its original position by a stiff spring. The strong spring forces act to hold the nodes in place and thus impose a no-slip boundary condition at the wall. The force exerted on wall node X w by the spring connecting it to its neighbor X w ′ is where d w is the initial spacing between the wall nodes. The force exerted on a wall node X w by the spring tethering it to its fixed initial position X 0w is The biofilm nodes at the bottom of the biofilm are initially placed along the wall co-linearly with the wall nodes. For the node X m located initially along the wall, a stiff tether force connects it to its initial position X 0m In the simulation of the biofilm in fluid flow, the flow is allowed to develop into a steady-state before allowing the biofilm to deform. This is accomplished by temporarily tethering each biofilm node X m to its initial position X 0m with a stiff force Once the flow has fully developed, the tethering force (Eq. 12) is released, thereby allowing the biofilm to deform. The spring constants appearing in Eqs. 9, 10, 11, and 12 are much stiffer than the springs connecting the biofilm nodes, i.e., A, B, C, D ≫ Ed with d a typical biofilm spring resting length.
The forces acting on the biofilm and wall nodes can be used to define a force denstity function on the biofilm-wall Lagrangian structure. With F * m the net viscoelastic force acting on biofilm node X m (from Eqs. 7, 11, and 12) and F * w the net spring force acting on wall node X w (from Eqs. 9, 10), the force density on the immersed structure is where the Dirac delta functions distribute the forces to the Lagrangian nodes at which they act. Inserting Eq. 13 into Eq. 3 gives the total body force acting on the fluid due to its interaction with the biofilm-wall structure as 26 which may be regularized on the discrete Eulerian grid via Eq. 5. In Supplementary Information-1, we discuss the origin of Eq. 13 and the calculation that leads to Eq. 14 via Eq. 3. We also show the connection between this formalism, which allows an a priori determination of the viscoelastic parameter values E and η, and other approaches that define a force density on the immersed structure through a "stiffness constant" with dimensions of force density. In Supplementary Information-2 and -3, we verify the accuracy of Eq. 14 by simulating a creep test on the biofilm structure where the elastic modulus E and dashpot viscosity η were derived directly from measured values of the shear modulus and relaxation time in biofilm experiments 10 . The creep test analysis indicates that the simulated elastic shear modulus and relaxation time are consistent with chosen experimental values.

Chemical Dynamics and Biofilm Growth Model
The spatiotemporal dynamics of the nutrient concentration in a fluid enivronment is governed by an advection-diffusion-reaction equation, where c is the concentration of a single nutrient expressed in Eulerian coordinates x and D c is the homgeneous diffusion constant. In the present model, the reaction term r accounts for the uptake of nutrient concentration defined in the Eulerian domain by the biofilm described as a Lagrangian immersed structure. Bacterial growth resulting from nutrient consumption can be modeled using Monod kinetics 47,48 , which defines a specific bacterial growth rate as a function of nutrient concentration, where µ max is the maximum specific growth rate, C is the local nutrient concentration expressed in the Lagrangian coordinates, and K is the half saturation constant. The time scale of our simulations is much shorter than the characteristic time of bacterial cell death 49 , thus a decay term is not included in Eq. 16. The net biofilm mass increases in time according to where M is the total biomass and is the specific growth rate averaged over the Lagrangian biofilm structure. The biomass accumulated during biofilm growth is distributed along the biofilm boundary according to a probability density function proportional to the local specific growth rate, where B is the portion of the Lagrangian biofilm boundary in contact with the surrounding fluid. As the biofilm grows it consumes nutrients within a region E r of the Eulerian domain that coincides with the Lagrangian biofilm structure. The reaction term in Eq. 15 is thus given by where dm dt is the net rate of change in nutrient mass due to consumption by the biofilm and is the local specific growth rate expressed in Eulerian coordinates. Finally, the net rate of biomass accumulation is proportional to the net rate of nutrient consumption, where 0 ≤ Y ≤ 1 is the yield coefficient. In the next subsection we describe our implementation of the Eqs. 15-22 within the IB2d software, where we modified the advection-diffusion upwind scheme to accommodate the reaction term.

Growth and Reaction Algorithm
The growth and reaction algorithm begins by initializing the Lagrangian biofilm structure and the nutrient concentration values defined in the Eulerian fluid domain. The biofilm is constructed by specifying the semicircular shape of the biofilm with radius r c and the distance between biofilm nodes ds. Then the software package DistMesh 50 is used to create a triangular mesh 25 that consists of Q = Q 0 nodes, N = N 0 triangles of average area a 0 , and S = S 0 connections (or edges). DistMesh uses the Delaunay triangulation algorithm to produce a mesh that contains mostly equilateral and uniform triangles. Existing nodes and the connections between them are maintained throughout the simulation, and additional nodes and connections are added during the biofilm growth as discussed below. The two-dimensional biofilm is taken to represent a χ = 1µm-thick slice of a three-dimensional biofilm. Thus, the bacterial biomass assigned to an average triangle is m 0 = λ v 0 n 0 χa 0 where λ is the mass density of a bacterium, v 0 is the volume of a bacterium, and n 0 is the number density of bacteria in a three-dimensional biofilm.
The mth node of the biofilm Lagrangian structure has the coordinate X m = (X m ,Y m ).
The concentration values are defined in the Eulerian fluid domain, which is discretized into a square grid with side length h. The grid point x i j = (x i , y j ) has the concentration value c i j . The concentration value C m at each Lagrangian biofilm node X m is determined by averaging the concentrations at the corners of the Eulerian grid square containing X m , where The concentration values on the Lagrangian biofilm structure can be used to compute the increase of biomass. The local specific growth rate, Eq. 16, is approximated as while the average specific growth rate, Eq. 18, becomes During a given time step of size δt, the incremental increase in biofilm mass δ M is found by solving Eq. 17 withΓ from Eq. 25 using the forward Euler method, where Nm 0 = M is the mass of the biofilm at the beginning of the time step. As the biofilm mass increases, the numbers of nodes and triangles increase. For each increase of biofilm mass δ M, the increase in the triangle number is δ N = δ M/m 0 . At the end of each time step, the change in triangle number accumulates from ∆N to ∆N + δ N. A number N ′ = ⌊∆N + δ N⌋ new triangles of area a 0 are then spawned at the boundary of the biofilm and the triangle accumulation is reset to the residue ∆N + δ N − N ′ . The growth of new triangles involves the addition of Q ′ new nodes. Each new triangle is added to the biofilm boundary using the probabilistic growth rule, Eq. 19, discretized as where p m is the probability the new triangle will be spawned at X m . A position X m ′ is chosen as one vertex of a candidate triangle (which is checked against conditions discussed below) by generating a random number ξ distributed uniformly on the interval (0, 1) and choosing an integer m ′ to satisfy Once the position X m ′ is selected, a second boundary point X m ′′ adjacent to X m ′ in the clockwise direction is selected (unless X m ′ is the last clockwise boundary node, then X m ′′ is adjacent in the counterclockwise direction). A third point X ′ is then computed as the point located along the bisector to the line joining X m ′ and X m ′′ such that the area of the isosceles triangle joining X m ′ , X m ′′ , and X ′ is a 0 . If X ′ is within a small distance ε from an existing biofilm boundary point X m ′′′ , then a new 6/20 triangle is formed by connecting X m ′ , X m ′′ , and X m ′′′ . If the candidate triangle overlaps the existing biofilm but the third point X ′ is not within ε of an existing point, or if the triangle overlaps the wall, a different position is chosen as the vertex of a new candidate triangle via Eq. 28. If neither of the former two conditions hold, a point X Q+n = X ′ , where n = 1, · · · , Q ′ , is spawned and connected to the two boundary points X m ′ and X m ′′ to form a new triangle. After each time step, the biofilm mass increases by an incremental amount δ M and it consumes an incremental mass of nutrient δ m = δ M Y . The nutrient is consumed within the region E r , which is computed as a polygonal region of the Eulerian domain that tightly encloses the biofilm Lagrangian structure using the MATLAB function 'polyshape'. The local specific growth rate defined in the Eulerian domain is discretized as The discrete form of the reaction term Eq. 20 is thus where h is the Eulerian grid size and δt is the time step. The growth and reaction algorithm can be summarized as follows: At the initial time t = t 0 , the biofilm structure is discretized as a triangular mesh containing Q = Q 0 nodes located at positions {X m } = {X 0m }, S = S 0 connections, and N = N 0 triangles with the triangle accumulation ∆N = 0 and the biomass assigned to an average triangle m 0 . The concentration field c i j is also initialized in the Eulerian domain at the coordinates x i j . To update the biofilm mass and reaction term at time t, we perform the following process:

Numerical Simulations
We simulate the deformation and growth of elastic and viscoelastic biofilms in fluid flow with a nutrient source. The biofilms are attached to the wall of a 2D microchannel, which has two walls parallel to the x-axis constructed from Lagrangian wall nodes. The length of the channel walls is L x and the distance between the walls is L y . A no-slip boundary condition on the walls is imposed by connecting adjacent nodes to each other and tethering each node to its initial position with stiff Stokes elements (Eqs. 9 and 10). Periodic boundary conditions are imposed at the inlet and outlet of the channel. The biofilm is initialized as a semicircular triangular-mesh structure located at the lower wall with Q 0 Lagrangian nodes, S 0 viscoelastic elements, and N 0 triangles. To anchor the biofilm to the channel wall, the biofilm nodes initially co-linear with the wall are tethered to their initial position with stiff Stokes elements (Eq. 11). The simulation set up is shown in Fig. 1(a).
We compute the deformation of each biofilm over a range of fluid flow rates. The fluid flows are simulated on a square-grid Eulerian fluid domain with the grid size h. Each flow is imposed by applying an external body force F ext = ρgx to each fluid node, where ρ is the fluid density, g is an acceleration constant, andx is the x-direction unit vector. The acceleration constant g is determined using the steady-state 2D Poiseuille flow solution such that a desired maximum fluid velocity u max occurs along the midchannel: To determine the accuracy of our flow, we first perform simulations in the microchannel without the biofilm structure by applying F ext to each fluid node with g determined from Eq. 31. The simulated maximum fluid velocity reaches its steady state after 20 seconds and matches the desired maximum fluid velocity u max (listed in Table 1) within 2%. Therefore, in simulations with the biofilm structure present, the biofilm is tethered in place for the first 20 seconds of the simulation by stiff Stokes elements (Eq. 12) as the flow evolves into a steady-state. Then the tethers are released allowing the biofilm to deform in response to the shear flow. We simulate the growth of each biofilm in different nutrient concentration configurations shown in Figs. 1(b)-(e) to account for the combined effect of the shear flow and nutrient source location on the growth of a viscoelastic or elastic biofilm. At 20 seconds the concentration is initialized as c 0 at the Eulerian grid points within the highlighted region of height h as indicated in Figs. 1(c)-(e). During the simulation, the concentration is held constant at c 0 on the same Eulerian grid points to represent a continuous supply of nutrient. The basecase of zero nutrient configuration Fig. 1(b) is used to test the deformation of the biofilm in the presence of flow but in the absence of growth. Then the growth and deformation of the biofilms are tested under the full-stream Fig. 1(c), downstream 1(d), and upstream 1(e) nutrient configurations. The numerical values of all simulation parameters are displayed in Table 1.

Elastic versus Viscoelastic Biofilms
We use our model to compare the behavior of a purely elastic biofilm composed of Stokes elements to a more biologically realistic viscoelastic biofilm composed of Maxwell elements in a shear flow. The elastic modulus is chosen to be E = 50 Pa in both cases. The viscoelastic biofilm has a dashpot viscosity η = 5 × 10 4 Pa.s, which is chosen to give a relaxation time near the 18 minutes that is common in biofilms over a wide range of elastic moduli 10 . In Supplementary Information-3, the creep test on the viscoelastic biofilm structure is presented to show how the elastic modulus, dashpot parameter, and relaxation time are related.
We first test the response of the elastic and viscoelastic biofilms to shear flow at various values of u max . A comparison 8/20 Bacterial mass density λ 1.12 g/ml 56 Volume of bacterium To quantify the extent of the biofilm deformation, we compute the x-and y-coordinates of the centers of mass, which are plotted versus time in Figs. 2(g) and (h). After a few hundred seconds, the x-coordinate (Fig. 2(g)) of the elastic biofilm changes little as the Stokes elements are nearly in equilibrium with the fluid flow, while the x-coordinate of the viscoelastic biofilm moves downstream at a nearly constant speed as the dashpots in the Maxwell elements continually increase in length. Relative to the x-coordinate, the y-coordinate (Fig. 2(h)) of both biofilms changes little over 2.5 hours. The change in the y-coordinate is about an order of magnitude smaller than the x-coordinate, which indicates the flow acts primarily to shear the biofilms.
Next, we test the growth of the elastic biofilm and the viscoelastic biofilm in the presence of the four nutrient configurations shown in Figs. 1(b)-(e) for values of u max = 0, 1, 2, 3, 4, 5 × 10 −6 m/s. The results of the simulation for the four concentration configurations under a flow value of u max = 5 × 10 −6 m/s at 2.5 hours are shown in Fig. 3. In each case, the viscoelastic biofilm grows further downstream than the elastic one. Zoomed-in snapshots of the biofilm shapes in Fig. 3 are shown in Supplementary Information-4.
The difference between the amount of growth of the elastic biofilm and the viscoelastic biofilm in our simulations is caused by the interplay between their growth toward the nutrient source and their downstream deformation in response to the shear   flow. To quantify the difference in the net growth between the two biofilms, we take the ratio of their biomass at 2.5 hours for each of the four nutrient configurations under each flow speed. The results, shown in Fig. 4(a), indicate the simulations with the elastic biofilm overestimate the growth when the nutrient source is located upstream and underestimate the growth when the nutrient source is located downstream. This is an intuitive result because the shear flow tends to push the viscoelastic biofilm into the vicinity of the downstream nutrient source more than the elastic biofilm (Figs. 3(c) and (g)), while the elastic biofilm resists moving downstream and remains in the vicinity of the upstream nutrient source more than the viscoelastic biofilm (Figs.

3(d) and (h)
). The greater growth of the elastic biofilm in the upstream configuration and the greater growth of the viscoelastic biofilm in the downstream configuration tend to cancel each other when the nutrient source is distributed symmetrically in the full-stream (Figs. 3(b) and (f)) and in the no-concentration case (Figs. 3(a) and (e)), giving a ratio near unity ( Fig. 4(a)). The growth and deformation of the biofilms results in the net migration of their centers of mass along the bottom of the channel. To measure the migration rate, we fit a least squares line to the x-coordinate of the center of mass versus time, shown in Fig. 4(b), and a least squares line to the y-coordinate of the center of mass versus time, shown in 4(c). In Supplementary  Information-4 we show the linear fits over 2.5 hours with the residuals displayed to verify their accuracy. The results for the x-direction migration rate ( Fig. 4(b)) indicate the viscoelastic biofilm always migrates further downsteam than the elastic biofilm. This result is intuitive because the viscoelastic biofilm always tends to deform more downstream in response to the fluid flow more than the elastic biofilm. In the case where the nutrient source is located upstream, the net migration of the biofilm can be in the upstream direction. At small flow speeds the viscoelastic biofilm migrates upstream, while at high flow speeds it migrates downstream. The transition from upstream to downstream migration in the upstream configuration does not occur for the elastic biofilm, as it migrates upstream at each flow rate. The migration rate of the two biofilms is also different in the full-stream nutrient configuration even though their net growth is approximately equal (Fig. 4(a)).
The fact that the x-direction migration rate downstream appears graphically as an approximately vertical shift from the x-direction migration rate upstream for both biofilms in 4(b) suggests that the migration rate due to growth and the migration rate due to shear deformation is approximately additive, and the migration rate due to growth is approximately constant. Assuming the upstream velocity is U up ≈ −U growth +U shear and the downstream velocity is U down ≈ U growth +U shear , where U growth and U shear are the magnitude of the contribution to the migration rate due to growth and the shear flow respectively, then U growth ≈ (1/2)(U down −U up ). Thus, taking the mean value of U growth over the range of u max for the viscoelastic and elastic biofilm gives a measure of for the migration rate due to growth: U growth = (0.021 ± 0.002) mm/hr where the reported error is the standard deviation of U growth . The migration speed due to the shear flow, U shear , can also be estimated from the data as U shear ≈ (1/2)(U down +U up ) ≈ U full , where U full is the migration speed of the biofilm in the full-stream nutrient configuration. The migration speed due to the shear flow at u max = 5 × 10 −6 m/s is U shear ≈ 0.056 mm/hr for the viscoelastic biofilm and U shear ≈ 0.016 mm/hr for the elastic biofilm. The y-coordinate migration rate, shown in Fig. 4(c), is similar for both biofilms at small flow rates, while at large flow rates, the elastic biofilm overestimates the y-direction migration rate.

Viscoelastic Biofilms
Measurements of biofilms reported in the literature indicate that over eight orders of magnitude of the elastic shear modulus (10 −2 − 10 6 Pa) biofilms have a common viscoelastic relaxation time of about 18 minutes 10 . We use our model to compare the growth and deformation of viscoelastic biofilms with relaxation times near 18 minutes over a range of two orders of magnitude in elastic modulus: 10 − 1000 Pa. The range was chosen to give reasonable numerical results for the values of shear flow used in our simulations, u max = 0, 1, 2, 3, 4, 5 × 10 −6 m/s over 2.5 hour-duration. For these flow rates, biofilms softer than 10 Pa have strains so large that numerical errors could become significant, while biofilms stiffer than 1000 Pa have strains so small that differences between biofilms became negligible. Moreover, large strains in our simulations are undesirable because biofilms are known to exhibit strain hardening at large strain values 41 , which our model does not capture. Nevertheless, the range of elastic moduli we simulate reliably capture a typical range of 70-700 Pa measured in microfluidic experiments 41 .
The relaxation time of the biofilm τ is related to the elastic modulus E and the axial viscosityη axial (i.e., the viscosity measured from axial stress versus axial strain) as In the Supplementary Information-2, we derive equations Eqs. 32 and 33 and discuss the quantity η water = 1.02 × 10 3 Pa.s, which is the drag coefficient of a fluid node in the simulation divided by the discretization length of the Lagrangian biofilm structure. In earlier studies 24,28 , the viscosity of the water alone was used to simulate the viscous response of biofilms without accounting for viscosity in the biofilm structure. We note that for a value of η = 5 × 10 4 Pa.s, η water accounts for a 10% correction toη axial . Thus, for values of E > 50 Pa, the viscosity of the water contributes less than 10% percent to the dashpot viscosity η determined by requiring τ = 18 minutes. We test the deformation of the viscoelastic biofilms in the response to shear flow without nutrient present ( Fig. 1(b)) for the range of elastic modulus values E = 10, 50, 100, 500 Pa, at a flow rate of u max = 5 × 10 −6 m/s. A comparison of the deformation of the biofilms is shown in Fig. 5. The results indicate that softer viscoelastic biofilms undergo larger strains and develop larger protrusions than stiffer biofilms (Figs. 5(a)-(d)). For each biofilm, the x-coordinate of the center of mass moves in the direction of the fluid flow (Fig. 5(e)). The closer to zero the asymptotic slope of the x-coordinate versus time becomes, the larger the elastic modulus is. It is expected that viscoelastic biofilms with large values of E behave like rigid biofilms with near zero strain rates ( Fig. 2(g)). Increasing E at a given shear stress reduces the strain of the biofilm (Eq. 7), and reducing the strain at fixed relaxation time (τ ∼ η/E for large E) reduces the strain rate of the biofilm (Eq. 8). The change in the y-coordinate of the center of mass over 2.5 hours is about an order of magnitude smaller then the x-coordinate of the center of mass for all the viscoelastic biofilms ( Fig. 5(f)), which indicates the flow acts primarily to shear the biofilms.
Next, we simulate the growth and deformation of the viscoelastic biofilms in the presence of the three nutrient configurations shown in Figs. 1(c)-(e) for values of E = 10, 50, 100, 500, 1000 Pa at u max = 5 × 10 −6 m/s. For other parameter values, see Table 1. Images of the results are shown in Fig. 6 for the full-stream, downstream, and upstream configurations. Figs. 6(a)-(d) suggest all viscolelastic biofilms grow a similar amount in the full-stream nutrient configuration. The softest biofilm appears to grow more in the downstream configuration 6(e) as compared to the upstream configuration 6(i), while the stiffest biofilm appears to grow more in the upstream configuration 6(l) as compared to the downstream configuration 6(h). Zoomed-in snapshots of the biofilms in Fig. 6 are shown in Supplementary Fig. 3.
To compare how the nutrient configuration affects the growth of the different viscoelastic biofilms in shear flow, we compute the biomass of each biofilm at 2.5 hours relative to the initial biomass for u max = 5 × 10 −6 m/s. The results, shown in Fig. 7(a), indicate that soft biofilms grow more in a downstream nutrient configuration while stiff biofilms grow more in an upstream configuration. Thus, there is a crossover, which occurs at E ≈ 80 Pa, shown as the vertical dashed line in Fig. 7(a). The crossover value of E is not universal but depends on the specific values of flow rate u max and the nutrient source concentration c 0 . The greater growth of soft bioflims in the downstream case and the greater growth of the stiff biofilms in the upstream case tend to cancel when the nutrient is symmetrically distributed as in the full-stream and no-concentration configurations, which leads to equal net growth over the full range of elastic moduli.
As seen in Fig. 6, when the viscoelastic biofilms grow, they migrate horizontally along the bottom of the channel and vertically away from the wall. To measure the horizontal and vertical migration rates, we fit least squares lines to the xcoordinate of the center of mass versus time and the y-coordinate of the center of mass versus time respectively for each biofilm in each nutrient configuration. In Fig. 7(b)   U down ≈ U growth +U shear , the upstream migration rate U up ≈ −U growth +U shear , and the full-stream migration rate U full ≈ U shear as discussed in the previous section. For all values of E, the difference between the downstream rate and the upstream rate gives a measure of the speed with which the biofilm migrates toward the nutrient, U growth ≈ (1/2)(U down −U up ) = (0.024 ± 0.002) mm/hr, which agrees with the result in the last section obtained by averaging U growth over the range of u max . At the elastic modulus value E = 10 Pa for the softest biofilm, U shear ≈ 0.112 mm/hr, and at E = 1000 Pa for the stiffest biofilm, U shear ≈ 0.009 mm/hr. In the limit of large values of E, the biofilm undergoes negligible strains due to the shear flow and thus U shear → 0 and the migration rate is dominated by U growth : in the upstream configuration U up → −U growth , in the downstream configuration U down → U growth , and in the symmetric full-stream configuration U full → 0. In Fig. 7(c), the vertical migration rate is plotted versus E in the downstream, full-stream, and upstream nutrient configurations. For small values of E, the biofilm grows in the downstream configuration similarly to the full-stream configuration, while for large values of E, the biofilm grows in the downstream configuration similarly to the upstream configuration. These behaviors can be understood as follows. For small values of E the biofilm is strained most downstream by the shear flow. In the downstream configuration, this large deformation means the biofilm is surrounded by the nutrient concentration, thus its local environment is like that of the biofilm in the full-stream configuration. In the upstream configuration, the biofilm is pushed out of the region of highest concentration, thus it grows more slowly. For large values of E, the biofilm is barely strained by the shear flow, thus the biofilm in the upstream configuration and the biofilm in the downstream configuration are surrounded by similar amounts of nutrient.

Discussion
In this work, we present a new model to simulate the growth of elastic and viscoelastic biofilms in shear flow under different nutrient source configurations. The interaction between the biofilm structure and the fluid is computed using the Immersed Boundary Method (IBM) and the growth is simulated using a probablistic rule based on local nutrient concentrations. Simulations are completed using the open source IB2d software [42][43][44] , which we modified to accommodate a reaction term in the advection-diffusion equation solver and to implement the growth scheme. The viscoelastic and chemical kinetic parameters in the simulations are taken directly from experimental values reported in the literature (Table 1).
Our method of applying IBM resolves the challenge 25 of determining viscoelastic parameter values a priori rather than having to perform a parametric search to match simulations to experiments. In IBM, the point forces exerted on the nodes of the Lagrangian structure must be transformed into a force density function. Typically this is accomplished (for an elastic structure) by defining "stiffness constants" that have dimensions of force per generalized Lagrangian volume. The difficulty is that these constants do not have a clear interpretation in terms of the actual spring constants or the elastic moduli. In a square lattice, elementary theory gives the interatomic spring constant k in terms of the elastic modulus E and the interactomic spacing d as k = Ed, independently of the dimensionality of the lattice 59 . Thus, it should be possible, given a Lagrangian structure with spacing d and an elastic modulus E, to a priori determine k. To accomplish this, we define the spring force between nodes using the spring constant k = Ed and use Dirac delta functions to define force densities directly from the point forces acting on the nodes of the Lagrangian structure. This procedure enables an exact evaluation of the force spreading equation (Eq. 3), which results in body forces in the Eulerian domain that take the same form as those employed in other computational fluid dynamics methods (Eq. 14), for example in the Method of Regularized Stokeslets 60 . The technical details of our approach are given in Supplementary Information-1 and -2, and the resulting simulation is validated in Supplementary Information-3.
Many earlier models treated biofilms as purely elastic structures, though biofilms are known to be viscoelastic. Thus, we use our model to find differences between simulations of a purely elastic biofilm and of a viscoelastic biofilm with the same elastic modulus of 50 Pa. In simulations with shear flow but without a nutrient source, the viscoelastic bifiolm develops a morphological structure with protrusions that the purely elastic biofilm does not develop. In simulations with shear flow and an asymmetric nutrient source, the elastic biofilm underestimates or overestimates growth rates depending on the configuration of the nutrient supply in the channel. With the nutrient supply located upstream from the biofilm, the elastic biofilm overestimates the growth rate. With the nutrient supply located downstream, the elastic biofilm underestimates the growth rate. In all nutrient configurations, even in the symmetric full-stream and no-concentration configurations, in which the growth rates are equal, the viscoelastic biofilm migrates more downstream than the elastic biofilm. In general, the differences between the elastic and viscoelastic biofilms are larger for higher flow rates. These results should be useful for future efforts to model biofilms using immersed structures because they illustrate some limitations of using a purely elastic structure.
Our simulations capture the biological reality that biofilms possess a common viscoelastic relaxation time of about 18 minutes over a wide range of elastic moduli 10 . We found that an interplay between biofilm deformation due to shear flow and biofilm growth toward a nutrient source resulted in a crossover value of the elastic modulus. Biofilms with elastic moduli below the crossover value grew faster when the nutrient source was located downstream than they did when the nutrient source was located upstream. Conversely, biofilms with elastic moduli above the crossover value grew faster when the nutrient source was located upstream than they did when the nutrient source was located downstream. We also found that biofilms migrated along the bottom of the channel wall. Soft biofilms migrated downstream in all nutrient configurations but migrated faster when the nutrient source was located downstream. Stiff biofilms migrated downstream when the nutrient source was located downstream but migrated upstream when the nutrient source was located upstream. Our results agree with the biophysical principle that by tuning their viscoelastic parameters, perhaps by modulating their polysaccharide production 7, 16, 18 , biofilms can adapt to their local nutrient conditions in shear flow to maximize their growth rate and their migration rate toward a nutrient source. In this work, we focus on the interplay between biofilm growth and deformation over a few hours. In the future, to further capture real characteristics and evolution of a viscoelastic biofilm, we will extend our model to include spatial porosity and heterogeneous structure and to exhibit detachment phenomena.