Incorporating Cellular Stochasticity in Solid–Fluid Mixture Biofilm Models

The dynamics of cellular aggregates is driven by the interplay of mechanochemical processes and cellular activity. Although deterministic models may capture mechanical features, local chemical fluctuations trigger random cell responses, which determine the overall evolution. Incorporating stochastic cellular behavior in macroscopic models of biological media is a challenging task. Herein, we propose hybrid models for bacterial biofilm growth, which couple a two phase solid/fluid mixture description of mechanical and chemical fields with a dynamic energy budget-based cellular automata treatment of bacterial activity. Thin film and plate approximations for the relevant interfaces allow us to obtain numerical solutions exhibiting behaviors observed in experiments, such as accelerated spread due to water intake from the environment, wrinkle formation, undulated contour development, and the appearance of inhomogeneous distributions of differentiated bacteria performing varied tasks.


Introduction
Bacterial biofilms provide basic model environments for analyzing the interaction between mechanical and cellular aspects of three-dimensional self-organization 1 arXiv:2401.07088v1[cond-mat.soft]13 Jan 2024 during development.Biofilms are formed when bacteria encase themselves in a hydrated layer of self-produced extracellular matrix (ECM) made of exopolymeric substances (EPS) [Flemming(2010)].This habitat confers them enhanced resistance to disinfectants, antibiotics, flows, and other mechanical or chemical agents [Hoiby (2010)].
Research on modeling biofilms has increased steadily during the past few decades resulting in the understanding of a number of features.Continuous models for uniform cell distributions are useful in basic culture systems [Stewart(2002)].Individual based models [Strorck(2014), Grant(2014)] and cellular automata [Laspidou(2004)] may capture variable thickness, density, and structure.However, current models focus more on deterministic mass transfer and extracellular structure, than in random cell processes.Interest on fluctuations in intracellular concentrations, for instance, has arisen due to their significance in phenotypic variability as well as in gene regulation and stochasticity of gene expression [kaern(2005), Wilkinson(2009)], with consequences for development and drug resistance [Birnir(2018)].
Recent experiments with Bacillus subtilus biofilms on agar provide a case study in which we can test models incorporating new aspects.Once bacteria adhere to a surface, they differentiate in response to local fluctuations created by growth, death, and division processes, to variations in the concentrations of nutrients, waste, and autoinducers, to cell-cell communication [Chai(2011)].Some of them become producers of exopolymeric substances (EPS) and form the extracellular matrix (ECM).EPS production increases the osmotic pressure in the biofilm, driving water from the agar substrate and accelerating spread [Seminara(2012)].In addition, the matrix confers the biofilm elastic properties.Wrinkles develop as the result of localized death in regions of high cell density and compression caused by division and growth [Asally(2012)].As the biofilm expands, complex wrinkled patterns develop, see Figure 1.This phenomenon is linked to gradients created by heterogeneous cellular activity and water migration [Espeso(2015)].Eventually, the wrinkles form a network of channels transporting water, nutrients, and waste to sustain it [Wilking(2013), Yan(2019)].Biofilm spread due to osmosis can be accounted for by two-phase flow models and thin film approximations [Seminara(2012)].Instead, wrinkle formation has been reproduced by means of Von Kármántype theories [Espeso(2015), Zhang(2016)].Delamination and folding processes are further analyzed in [benamar(2014)] by means of neo-Hookean models.In [Carpio(2018)], a poroelastic approach provides a unified description of liquid transport and elastic deformations in the biofilm.To incorporate fluctuations in a more natural way, here we propose a mixture model allowing to distinguish the different phenotypes forming the film.
Biofilm structure is greatly influenced by environmental conditions.When they grow in flows, we find bacteria immersed in large lumps of polymer, typically forming fingers and streamers [Drescher(2013)] in the surrounding current.In contrast, biofilms spreading on air-agar interfaces contain small volume fractions of extracellular matrix [Seminara(2012)], producing wrinkled shapes with internal water flow.This motivates different treatments of the extracellular ma- trix, see [Kreft(2001), Jayathilake(2017)] for biofilms in flows and [Strorck(2014), Grant(2014), Seminara(2012), Lanir(1987)] for biofilms on interfaces with air or tissues, for instance.In the latter case, when internal fluid flow is taken into account, the small fraction of matrix is usually merged in one biomass phase with the cells [Seminara(2012), Lanir(1987)].Some experimental studies suggest a viscoelastic rheology for biofilms [Shaw(2004), Charlton(2019)].The analysis of the mixture and poroelastic models we consider shows that, depending on the volume fractions of solid biomass and fluid, the viscosity of the fluid, the Lamé constants of the solid, the densities, and hydraulic permeability of the fluid/solid system, the characteristic time for variations in the displacement of the solid, and the characteristic length of the network in the macroscopic scale, the resulting mixture can be considered as monophasic elastic, monophasic viscoelastic, or truly biphasic mixture/poroelastic [Burridge(1981), Kapellos(2010)].
The paper is organized as follows.Section 2 introduces the solid-fluid mixture model.Section 3 discusses ways to incorporate details of cell behavior.We present a cellular automata approach based on dynamic energy budget descriptions of bacterial metabolism.With the aid of asymptotic analysis [Seminara(2012), Witeslki(1997)], we construct numerical solutions displaying behaviors consistent with experimental observations.Finally, Section 4 discusses our results, the advantages, and limitations of our approach, as well as future perspectives and possible improvements.

Mass Balance
Under the equipresence hypothesis of mixtures, each point x in the biofilm can be occupied simultaneously by both phases.In addition, we assume that no air bubbles or voids form inside the biofilm.If ϕ s (x, t) denotes the volume fraction of solid and ϕ f (x, t) the volume fraction of fluid, then ϕ s + ϕ f = 1.The standard densities of biological tissues and agar are similar to the density of water ρ f = 10 3 kg/m 3 (typical relative differences of order 10 −2 ).Thus, we will take the densities of all constituents and the mixture to be constant and equal to that of water [Seminara(2012)]: Then, the balance laws for the fractions of solid biomass ϕ s and fluid ϕ f are [Lanir(1987), Seminara(2012)] where r s (ϕ s , c n ) represents biomass production due to nutrient consumption, whereas v s and v f stand for the velocities of the solid and the fluid components, respectively.Biomass production can be accounted for through a Monod law , where c n is the nutrient concentration, K n is a constant that marks the onset of starvation, and k s the uptake rate, which can be approximated by 1+αm τ , τ being the doubling time for the specific bacteria and α m a factor representing polymeric matrix production [Seminara(2012)].
Adding up Equation (1), we find a conservation law for the growing mixture: where v = ϕ s v s + ϕ f v f is the composite velocity of the mixture and is the filtration flux, w being the relative velocity.

Driving Forces
Forces inducing motion are of a different nature: inner stresses, inertial forces, interactive forces, and external body forces.We discuss here the constitutive relations and fluxes for incompressible solid-fluid mixtures in the case of infinitesimal deformations and under isothermal conditions [Lanir(1987)].

Stresses in the Solid and the Fluid
When a large number of small pores are present in the biofilm, the stresses in the fluid are where p is the pore hydrostatic pressure.In the presence of large regions filled with fluid, the overall fluid shear stresses should be considered too.The stresses in the solid arise from the strain within the solid and from the interaction with the fluid.Under small deformations, and for an isotropic solid, we have where u s is the displacement vector of the solid biomass; ε(u) the deformation tensor; and λ s , µ s , the Lamé constants, related to the Young E and Poisson ν moduli by λ =

Interaction and Inertial Forces
In most biological samples, the velocities v f and v s are small enough for inertial forces to be negligible: ρ s a s ≈ ρ f a f ≈ ρa ≈ 0, where a s , a f , a represent the solid, fluid, and composite accelerations.Thus, we will work in a quasi-static deformation regime.
The interaction forces act on the two components.They are opposite in sign and equal in magnitude, as a result, their combined effect vanishes on the tissue.We consider two kinds: filtration resistance and concentration gradients in chemical potentials.
The filtration resistance arises from the interaction between fluid and solid particles.Per unit volume, these forces are ϕ s f s , ϕ f f f and satisfy ϕ s f s +ϕ f f f = 0.In the absence of inertial effects and concentration-viscous couplings f f = −αq, where α(ϕ f ) is the resistivity matrix and q the filtration flux.For isotropic elastic solids with isotropic permeability where k h = k µ f is the hydraulic permeability, k being the permeability of the solid, and µ f the fluid viscosity.Typically, , where ζ is a friction parameter often taken to be ζ = µ f ξ(ϕs) 2 > 0 and ξ represents the "mesh size" of the matrix network.
The concentration forces in the fluid ∇π f are , where Vf is the molar volume of the fluid and µ f,c is the concentration contribution to the chemical potential of the fluid µ f .Under isothermal conditions Similar relations hold for concentration forces ∇π s in the solid, which satisfy ϕ s ∇π s + ϕ f ∇π f = 0.

Equations of Motion
The theory of mixtures hypothesizes that the motion of each constituent is governed by the usual balance equations, as if it was isolated from the other one.Neglecting inertial terms, and in the absence of external body forces, the momentum balance for the solid and the fluid reads Using the expressions for the stress tensors (4) and ( 5), these equations become Combining ( 8), (6), and (3) we obtain the Darcy law in the presence of concentration gradients Adding up equations (8), we find an equation relating solid displacements and pressure div σs (u s ) − ∇p = 0. (10) The solid velocity is then v s = ∂us ∂t .These equations are complemented by the conservation of mass (1) and ( 2), which now reads as The flux (9) can be rewritten as q = k h div(− σs + π f I), where − σs + π f I is the swelling stress.In biphasic swelling theory [Wilson(2005)], it is customary to work with p − π f = p f , where p f is associated to the fluid chemical potential by ( 7) and π f is identified with the osmotic pressure created by the concentration of a specific chemical [Wilson(2005), Ghassemi(2003)].The osmotic pressure in the biofilm is caused by the concentration of EPS produced by the cells and can be taken to be proportional to the volume fraction of solid biomass π f = Πϕ s , Π being the osmotic compressibility [Seminara(2012)].Equation (10) then motivates the introduction of effective constitutive laws for the whole mixture of the form [Wilson(2005)] σ(u) = σs (u) − pI = σs (u) − (p f + π f )I, as usual in poroelastic theory.

Final Equations
Summarizing the main governing equations, we get in the region occupied by the biofilm Ω b (t), which varies with time.In equilibrium, q = v s = f f = 0.At the biofilm boundary, the jumps in the total stress vector and the chemical potential vanish: when applicable.In this quasi-static framework, the displacements u s depend on time through the motion of the biofilm boundary.
If we need to track the variations of the nutrient concentration, we may use effective continuity equations for chemical concentration in tissues [Chen(2007), Sacco(2107)].For the limiting concentration c n : where d n is an effective diffusivity [Wood(2002)].Setting d n,s and d n,f , the diffusivities in the biomass and liquid d n,f .The source r n (ϕ s , c n ) represents consumption by the biofilm, k n being the uptake rate and K n the half-saturation constant.Zero-flux boundary conditions are imposed at the air-biofilm interface.Instead, at the agar-biofilm interface, we may impose a constant concentration through a Dirichlet boundary condition.Being more realistic, we couple this diffusion equation to another one defined in the agar substratum Ω a (t) with zero source and transmission conditions at the interface [Espeso( 2015 )].
Solving Equations ( 12)-( 15), studying the evolution of a biofilm also requires tracking the dynamics of the biofilm boundary, see Figure 3.In principle, there are two boundaries of a different nature: the air-biofilm interface and the agarbiofilm interface.

Motion of the Air-Biofilm Interface
During the first stages of biofilm spread, the agar-biofilm interface remains flat, whereas the biofilm reaches a height x 3 = h(x 1 , x 2 , t), see Figure 3a.Integrating (2) in the x 3 direction we obtain ∂x3 dx 3 = 0, x1 , x2 and x3 being the unit vectors in the Cartesian coordinate directions.By Leibniz's rule . Inserting this identity in (17), we find the equation where To obtain a closed equation for the height h we need to calculate the velocity of the solid v s = dus dt , the modified pressure p f and the volume fraction of fluid from ( 12) and ( 13).This equation is able to describe accelerated spread due to osmosis, at least in simplified geometries, as we illustrate next.
From Equation ( 18), we derive an approximated equation for the early evolution of the height of a circular biofilm, see Appendix 4 for details and assumptions A simplified version has self-similar solutions.Restoring dimensions, they take the form Replacing (1 + 3/2) by 1 in (21), we recover the self-similar solution found in [Seminara(2012)], with gµ f instead of µ f (µ f being the fluid viscosity) and the Lamé coefficient of the solid biomass µ s instead of the viscosity of the fluid biomass µ s .
Figure 4 compares the time evolution of the biofilm height profiles starting from a smoothed version of (21).Notice that (21) only makes sense when R 2 > 3/2 r 2 , and that the slope diverges at r = 2/3R.Experiments show that a thin biofilm layer precedes the advance of the biofilm bulk [Seminara(2012)].We set h = h ∞ > 0 beyond that point.The dashed green line in Figure 4 represents the numerical solution of (20), with R given by ( 21) for K = 10 −5 , and h inf = 10 −3 , replacing (1 + 3/2) by 1 as in [Seminara(2012)].The dotted red line and the solid blue line depict the numerical solution of ( 20) and ( 19), respectively, with R given by ( 21) and keeping the same data.They all show the transition from vertical growth to horizontal spreading as time goes on.The effect of the additional time derivatives in ( 19) is to flatten the profiles.
When the interface biofilm/agar is not flat, but admits a parametrization of the form x 3 = ξ(x 1 , x 2 , t), as in Figure 3b,  20) and ( 19), respectively, with R given by ( 21) and keeping the same data.We can observe the transition from an initial stage in which increase in biofilm height dominates to a stage with faster horizontal spread.The green line is a reference self-similar approximation.
Repeating the previous computations in the interval [ξ, h], the equation for the biofilm height becomes Knowing ξ, this equation can be solved numerically coupled to ( 12) and ( 13).

Motion of the Agar/Biofilm Interface
Equations for the dynamics of the agar-biofilm interface follow using a Von Karman-type approximation, as the thickness of the biofilms is small compared to its radius.Although initially flat, the displacements in the direction orthogonal to the interface may become large.Thus, the linear definition of the strain and stress tensors in ( 5) is replaced by [Landau(2010)] which includes nonlinear terms, as well as residual strains ε 0 i,j .We denote the in-plane displacements by u = (u 1 (x 1 , x 2 , t), u 2 (x 1 , x 2 , t)) and the out-ofplane displacements of the interface by ξ(x 1 , x 2 , t).The coordinates (x 1 , x 2 ) vary along the 2D projection of the 3D biofilm structure on the biofilm/agar interface.Equation ( 14) becomes div σ = ∇(p f + π f ), with σ given by ( 5) and ( 23).Formally, this allows us to identify the biofilm with an elastic film growing on a viscoelastic agar substratum.The pressure terms become residual stresses.Then, the interface motion is governed by the equations [Huang(2006), Espeso(2015)]: where h a is the thickness of the viscoelastic agar substratum and µ v , ν a , and η a its rubbery modulus, Poisson ratio, and viscosity, respectively.The tensor σ is given by with ε defined in (23); ν and E being the Poisson and Young moduli of the biofilm (5), respectively; and σ 0 represents the residual stresses.The bending stiffness is D = Eh 3 12(1−ν 2 ) , h being the initial biofilm thickness.Here, the first Equation ( 24) describes out-of-plane bending ξ, and the second one (25) governs in-plane stretching for the displacements u = (u 1 , u 2 ).Modified equations taking into account possible spatial variations in the elastic moduli are given in [Iakunin(2018)].
To identify the relevant scales governing the evolution of the agar-biofilm interface we nondimensionalize ( 24) and (25).Making the change of variables where R is the approximate biofilm radius, and setting R = γh, the dimensionless equations become where Wrinkled structures develop when the nonlinear terms are large enough, therefore γ = R h must be large enough.
The residual stresses σ 0 in (26) can be estimated averaging the osmotic and fluid pressure contributions to the three-dimensional biofilm.If the solution of ( 12)-( 15) in the biofilm Ω b (t) is known, σ 0 ij could be estimated from where x 3 = h and x 3 = ξ define the two biofilm interfaces with air and agar, see Figure 3. Analytical approximations ( 19)-( 21) of the the biofilm height h in early stages of the biofilm evolution allow for simple simulations of the onset of wrinkle formation.Figure 5 uses these asymptotic profiles to compute pressures and velocities by means of ( 49)-( 52).Starting from an initially flat biofilm, ( 29) suggests that we should consider stress profiles of the form , which are nondimensionalized dividing by E. The first two terms reflect the stresses due to growing height and radius, whereas the last one accounts for the osmotic pressure.Inserting these residual stresses in ( 24)-( 26) we generate small inhomogeneities and wrinkles in Figure 5.However, these approximations neglect spatial variations in concentrations, as well as changes in cell behavior, and therefore, in stresses and pressures.Therefore, the patterns display soon an unrealistic behavior, with wrinkles excessively growing in the central region.
Solving the full set of coupled equations we have derived is very costly and faces severe numerical difficulties at contact points.Alternatively, we may set σ 0 = 0 and work with the residual strains ε 0 in ( 23), which can be related to growth tensors created by stochastic cell processes as we discuss next.

Incorporating Cellular Behavior
Cells within a biofilm differentiate to perform different tasks, and can deactivate due to lack of resources or die as a result of biochemical stress and waste accumulation [Asally(2012), Chai(2011)].Such variations in the biofilm microstructure affect the overall shape [Espeso(2015)].Cell activity enters the previous deterministic model through the biomass creation term g(c n )ϕ s in (12), the nutrient consumption term r n (ϕ s , c n ) in ( 16), and the residual stresses σ 0 in (26).However, this does not account for cell death, cell deactivation and cell differentiation.
Differentiation implies changes in phenotype while preserving the same genotype.For B. subtilis biofilms, the differentiation chain through which different cell types originate is established in [Chai(2011)], see Figure 6.Initially, we have a population of similar alive cells glued together in a matrix, most of which have lost their individual motility.All of them secrete ComX.If the concentration of ComX becomes large enough, some cells differentiate and start producing surfactin, losing their ability to reproduce.For large enough surfactin concentrations, other normal cells differentiate and become EPS producers.These cells reproduce more slowly than normal cells.Cells may also die due to biochemical stresses [Asally(2012)], preferentially at high-density regions, such as the agarbiofilm interface.In the upper regions of the biofilm, depletion of resources may trigger deactivation of cells, which become spores.Undifferentiated cells retaining some motility are restricted to the bottom edges [Chai(2011)].
To a large extent, these processes have a random character.Hybrid models combine stochastic descriptions of cellular processes with continuous equations for other relevant fields.This allows us to consider the inherent random- ness of individual bacterial behaviors as well as local variations [Chai(2011), Mehta(2008)].We will explain how to introduce cell variability in the mixture model next.

Cellular Automata and Dynamic Energy Budget
In a cellular automata approach, space is divided in a grid of cubic tiles.This grid is used to discretize all the continuous equations: concentrations, deformations, pressures, etc.To simplify geometrical considerations, we initially assume that each tile can contain one bacterium at most.We describe bacterial metabolism using a dynamic energy budget framework [Birnir(2018), Kooijman(2008)].According to this theory [Kooijman(2008)], cells create energy from nutrients/oxygen, which they use for growth, maintenance, and product synthesis.Damage-inducing compounds can cause death.The metabolism of each cell C j is described by the system: dve,j dt = (1 − α)r e,j v j , r e,j = kr j + k ′ , dqj dt = e j (s G ρv j q j +h a )(ν ′ −r j ) − (r j +r e,j )q j , dhj dt = q j − (r j + r e,j )h j , where the cell variables are their energy e j (t), volume v j (t), produced volumes of EPS v e,j (t), acclimation a j (t), damage q j (t), hazard h j (t), and survival probability p j (t), for j = 1, . . ., N , N being the total number of cells.These equations are informed by the value of continuous fields at the cell location x (the tile of the cellular automata grid containing the cell): nutrient concentration c n (x, t), polymeric matrix concentration c e (x, t), surfactin concentration c s (x, t), ComX concentration c cx (x, t), and environmental degradation ε(x, t), which are gov-erned by The parameters ν, m, g, a M , and ρ are the energy conductance, the maintenance rate, the investment ratio, the target acclimation energy, and mass density for the bacteria under consideration, respectively.Other coefficients are the multiplicative stress coefficient s G , the maintenance respiratory coefficient ν m , the noncompetitive inhibition coefficient K v , and the environmental degradation coefficients γ and ν ε .The parameters d n , d cx , d s , d e , and d ε stand for diffusivities.The Dirac masses δ j are equal to 1 at the location of cell C j and zero outside.
The production rates r s,j and r e,j are zero, except when the cell is a surfactin producer, or an EPS producer, respectively.In the latter case, r e,j = kr j + k ′ , where k and k ′ correspond to constants controlling the chemical balances for polymer production.The parameter α ∈ [0, 1] regulates the fraction of produced polymer that remains in a monomeric state and diffuses as c e , instead of becoming part of larger chains that remain attached to the cells forming the matrix v e .
In this framework, bacteria C j die with probability 1 − p j .Taking a cellular automata view, we modify the cell nature according to selected probabilities, which are defined in terms of concentration values at the cell location.A normal bacterium becomes a surfactin producer with probability p cx = The simulation started from a circular seed with a diameter of 60 cells, and nonuniform height.Each colored box in the slices represents one cell.The brown areas representing dead cells appear at the bottom of three initial peaks.We set knL 2 dnKn = 0.01, kcxL 2 dcxKcx = 0.01, and ksL 2 dsKs = 0.8., where L is a reference length representing the tile size (approximately the bacterium size 2 µm).
In principle, when a bacterium divides, the daughters occupy the space left by it, while pushing the other bacteria.Dealing with the geometrical aspects of arrangements of dividing bacteria is a complicate issue for which different approaches have been explored [Strorck(2014), Grant(2014)]; it is out of the scope of the present work.For simplicity, we consider here that space is partitioned in a grid of cubic tiles, as explained earlier, and this grid is used to discretize all the continuous equations for concentrations, displacements, pressures, etc.Each tile may contain at maximum one bacterium, the tile size is the maximum size a bacterium may attain.Once a bacterium divides, one daughter remains in the tile, whereas the other occupies a random neighboring tile, either empty, or containing a dead cell, which is reabsorbed.In the absence of them, it will push neighboring cells in the direction of minimal mechanical resistance, that is, minimal distance to air.The resulting collection of tiles defines the new Ω b .Figure 7 represents the evolution of a growing biofilm seed considering only division processes with c n (0) = K n and knL 2 dnKn = 8.Notice the development of contour undulations.The resulting full computational model would proceed as follows.
• Each bacterium is initially classified as normal, surfactin producer, EPS producer, or inert.Bacteria are distributed in the tiles x of the grid.The empty space around them is filled with water and dissolved substances.
In this way, we may compute the volume fractions of biomass ϕ s (x, 0) and fluid ϕ f (x, 0) in each tile x, as well as the osmotic pressure Π(x, 0).The pressure p(x, 0) is obtained from (13) with v s = 0 and σ 0 (x, 0) from ( 29).
• We compute stationary solutions of the Equation ( 31) for v f = 0 by a relaxation numerical scheme.All except the equation for c n are solved using the grid defining Ω b (0) with no flux boundary conditions.The equation for c n is solved in the biofilm-agar domain, that is, Ω b (0) ∪ Ω a (0), imposing continuity of concentrations and fluxes at the agar-biofilm interface and no flux boundary conditions at the air-biofilm interface.
Evolution with a time step dt: From time t − dt to t.
• We update ξ(t) using ( 24)-( 29).Keeping the grid tile fixed, we shift the contains of the tiles upwards of downwards to reflect the evolution of ξ(t) when the deflections are large enough.
• We update the cellular fields solving (30) with the stationary concentration values for (31).The time derivatives in (31) are small, so that time evolution is driven by the source terms reflecting cell activity and changes in the biofilm boundaries.
• We update the bacterial status checking whether normal bacteria become surfactin or EPS producers, whether any of them deactivates or dies, and whether they divide, with the probabilities assigned to each situation.In case a bacterium divides, we reallocate the newborn cell.
• In the resulting biofilm configuration Ω b (t), we compute the volume fractions of biomass ϕ s (x, t) and fluid ϕ f (x, t) in each tile.This also provides the osmotic pressure π f (x, t).The fluid pressure p(x, t) is obtained from (13), the displacements u s (x, t) from ( 13), and σ 0 (x, t) from ( 29).The solid velocities are approximated by v s (x, t) ∼ (u s (x, t)−u s (x, t−dt))/dt.Then, the fluid velocity is given by (9).
• We yet need to take into account water absorption from agar.To do so, we solve ∂ϕ f ∂t + div(v f ϕ f ) = 0 in the biofilm/agar system.Alternatively, we can solve only in the biofilm, using ϕ f = 0 at the biofilm/agar interface and ∂ϕ f ∂n = h R ϕ f for boundary conditions, where h and R are reference values for the biofilm height and radius.Then, we revise the biofilm configuration, creating water tiles with probability ϕ f and shifting the contains of the neighbouring tiles.This provides the final biofilm configuration Ω b (t), that is, the occupied tiles, their contents, the bacterial status and fields, as well as the values of the continuous fields at each tile.This process mixes the stochastic evolution of some cellular processes with continuous equations for a number of fields.In case, any of the fields computed from the tile configuration is not smooth enough to solve the required partial differential equations, we filter it using a total variation based filter introduced in [Carpio(2018)] to avoid such artifacts.Figure 8 depicts the formation, coarsening, and branching of wrinkles in an spreading biofilm when the residual stresses are fitted by an empirical circular front approximation of magnitude −0.1 advancing one grid box every 14/τ seconds.
This hybrid model also introduces a number of parameters that should be calibrated to experimental data, not yet available.The parameters appearing in the dynamic energy budget equations have been fitted to experimental measurements for Pseudomonas aeruginosa biofilms under the action of antibiotics [Birnir(2018)]; fitting to Bacillus subtilis would require new specific experiments.

Balance Equation Approach
The macroscopic effect of the presence of differentiated bacteria can partially be understood by means of additional balance equations, inserting specific information in them.Let us set ϕ a and ϕ d as the volume fraction of active and dead cells, respectively.We introduce an additional volume fraction of inert cells ϕ i , in such a way that ϕ s = ϕ a + ϕ i + ϕ d .The balance equations become where r s (ϕ a , ϕ d , c n ) = g(c n )ϕ a − k r ϕ d , k r being the rate of reabsorption of dead cells and c w the concentration of waste.The concentration of nutrients still obeys ( 16), replacing ϕ s by ϕ a in the consumption term, whereas the concentration of waste c w obeys a similar reaction-diffusion equation with source r w (ϕ a ) = k w ϕ a , k w > 0. Here, g i (c n ) is positive for small enough values of c n and negative otherwise.For instance, we might take g i (c n ) = α − cn cn+Kn , α ∈ (0, 1).We assume that dead and alive cells move with the velocity of the solid biomass v s .Adding up Equations ( 32)-( 35), we recover the relations (2) and ( 13).The displacements of the solid biomass u s still obey ( 14) with two modifications.First, the osmotic pressure π f depends only on the fraction of cells producing EPS, which must be alive.Thus, π f = π(ϕ a ).Second, the elastic constants λ s and µ s may vary spatially in case necrotic regions containing a large density of dead cells or swollen regions appear.We focus here on the effect of necrotic regions on liquid transport within the biofilm.Figure 9 illustrates water accumulation in regions with an initially high volume fraction of dead cells.
To investigate the spatial distribution of cells secreting autoinducers, we introduce additional volume fractions of active cells ϕ a = ϕ u + ϕ surf + ϕ eps , where ϕ surf and ϕ eps stand for cells producing surfactin and EPS respectively, whereas ϕ u are undifferentiated active cells.The balance equations governing the different subpopulations are where If we consider dead and inert cells too, systems ( 32)-( 35) should be updated replacing g(c n )ϕ a by g(c n )ϕ u + g e (c n )ϕ eps in Equation ( 32) and in the definition of r s for Equation (35).Likewise, systems ( 36)-( 38) should be updated including transfer to and from the inert status.

Discussion
Growth of cellular aggregates involves mechanical, chemical, and cellular processes acting in different time scales.Bacterial biofilms provide basic environments to test hypotheses and mathematical models against experimental observations.Recent experimental work with Bacillus subtilis reveals a host of phenomena during biofilm formation and spread.Different approaches have been exploited to account for different aspects: thin film equations and two-phase flow models for accelerated spread caused by osmosis [Seminara(2012)], elasticity theory for the onset of wrinkle formation [Asally(2012), Iakunin(2018)], Von Karman-type approximations for wrinkle branching [Espeso(2015), Zhang(2016)], and Neo Hookean models for contour undulations [benamar(2014)] and fold formation.In principle, poroelastic models allow to consider liquid transport and elastic deformation in a unified way [Carpio(2018)], though detachment and blister formation require further developments [Yan(2019)].Current models take mainly a deterministic point of view, thus, random cell behavior linked to fluctuations is poorly accounted for.However, cell differentiation [Chai(2011)] to incorporate new phenotypes performing new tasks, such as autoinducer and EPS matrix production, plays a key role in biofilm development.Elementary cellular automata approaches were implemented in [Espeso(2015), Carpio(2018)] and used to generate nonuniform residual stresses partially defining the biofilm shape.Here, we develop a hybrid computational model, combining a solidfluid mixture description of mechanical and chemical processes with a dynamic energy budget based cellular automata approach to cell metabolism.
Cellular automata representations are convenient from a computational point of view, as they allow for simple rules to transfer information between individual cells and the film.However, they provide too crude a representation of bacterial geometry.In our framework, this representation could be improved by resorting to different agent based models.Individual-based models, originally developed to study biofilms in flows [Kreft(2001), Jayathilake(2017)], have recently been adapted to describe biofilms spreading over air-agar interfaces and solidsemisolid interfaces [Strorck(2014), Grant(2014)].Similarly, immersed boundary methods introduced to study bodies immersed in fluids are being extended to study biofilm spread in flows [Stotsky( 2016)] and at interfaces [Dillon(2008)].We could resort to Individual based or Immersed boundary approaches for a better description of bacterial geometry and their spatial arrangements.
Working with biofilms spreading on an air-agar interface, we have chosen to represent the presence of small fractions of polymeric matrix in an effective way, as done in previous related work [Seminara(2012), Lanir(1987)].The biomass formed by bacteria and polymeric threads is considered one phase [Seminara(2012)], with elastic properties as in [Lanir(1987)].The liquid transporting dissolved chemicals is considered a fluid phase.Production of EPS also affects internal liquid flow by osmosis, mechanism we include in our equations for the fluid phase.Depending on the relative fractions and the properties of each phase as well as the characteristic times and lengths, the whole system may display an elastic, fluid, viscoelastic, or truly poroelastic behavior [Burridge(1981), Kapellos(2010)].This formulation allows to derive effective equations for the dynamics of the interfaces including the effect of biomass growth, fluid, and osmotic pressures through residual strains and stresses.Resorting to individual-based or immersed boundary representations of cells, we might describe the polymeric matrix as a network of threads instead [Strorck(2014), Stotsky(2016)], but we should define heuristic rules for their behavior.
Constructing numerical solutions of the full model is a computational challenge, out of the scope of the present work.Instead, we construct numerical solutions, in particular, geometries, guided often by asymptotic simplifications.In this way, we show that the model is able to reproduce behaviors experimentally observed: accelerated spread due to water intake [Seminara(2012), Yan(2019)], wrinkle formation and branching [Asally(2012), Wilking(2013), Yan(2019)], layered distributions of differentiated cells [Chai(2011)], development of undulations in the contour [benamar(2014), Yan(2019)], and appearance of regions containing a high volume fraction of water [Wilking(2013), Yan(2019)].Existing models are devised to explain specific behaviors in relation with particular experiments.An advantage of our approach is that a single model can be used to display all those behaviors and to simulate or even analyze under which conditions they are observed, as the model allows for asymptotic analysis in specific situations.The partial study of different phenomena also suggests empirical expressions for magnitudes representing cellular activities required by the mixture model, such as source terms or residual stresses, which can be inserted in it to reduce computational costs.Our simulations of biofilm spread and wrinkle formation use parameter values experimentally measured for Bacilus subtilis biofilms in [Seminara(2012), Asally(2012)], producing reasonable qualitative and quantitative results.However, the parameters for the dynamic energy budget systems for cell metabolism, as well as those appearing in the concentration equations are taken from Pseudomonas aeruginosa studies [Birnir(2018)].The probability laws for the cellular automata model and the balance equations for differentiated cell populations involve additional unknown parameters.Thus, our model involves a collection of parameters that should be fitted to experimental data, specially as far as cell metabolism is concerned.Experimental measurements of bacterial dynamics allowing to fit such parameters are yet missing.by matching a two phase flow model for the biofilm with another one for agar, we keep it here.The solution consistent with these boundary conditions is , where the contribution ε −1 = R h disappears thanks to the boundary condition.For the expansion of ϕ in powers of ε p to be consistent we need ε p /ε << 1.Notice that u (0) 3 = x 3 t implies u (0) 3,x3 x3=1 = t.To compute u 1 and v 1 we use Equation ( 47), and the derivatives with respect to x 1 of p and ϕ, which enter this equation through the dependence on h(x 1 , t).To address this issue properly, we switch back to dimensional variables using P = gh 2 ζ and x ′ 3 = x 3 /h to get discarding in p the lower order terms which do not contribute to the derivatives.Then, Equation (47) yields Integrating twice and applying the boundary conditions (40) at 0 and h we find the displacement Differentiating the displacements u 1 and u 3 with respect to time we get the velocities Once we approximate expressions for velocities, pressures, and volume fractions are available in the x 1 x 3 plane, the equation for the free boundary (18) becomes h t + ∂ ∂x1 h ∇p, where v is the volume averaged velocity and v s the velocity of the solid, given by (52).Setting ϕ ∼ ϕ ∞ , we have v • x3 0 = gh (1−ϕ) 2 (1−ϕ∞) 2 ∼ gh, and the flux becomes Performing the change of variables h = e th , h t = e th t + e th , and dropping the symbol˜for simplicity, we find the desired equation.
In radial coordinates the same arguments work.The equation for the height is then: In dimensionless variables r = R 0 r, t = g −1 t, h = h 0 h.Dropping the symbol ˜for simplicity the equation becomes Performing the change of variables h = e th , h t = e th t + e th , and dropping again the symbol˜for simplicity, we find ( 19).

Figure 1 :
Figure 1: Virtual visualization of a biofilm spreading on agar.

2Figure 2 :
Figure 2: Schematic structure of a biofilm: (a) View of the macroscopic configuration: a biofilm on an agar-air interface.(b) Microstructure formed by biomass (polymeric mesh and cells) and fluid containing dissolved substances (nutrients, waste, and autoinducers).

Figure 3 :
Figure 3: Schematic representation of a biofilm slice, with moving air-biofilm interface h and agar-biofilm interface ξ.(a) Initial stages: ξ is flat.(b) Later evolution: ξ deviates out of a plane.

Figure 4 :
Figure 4: Biofilm height at dimensionless times 0 (a), 1 (b), 6 (c), and 7 (d) for K = 10 −5 and h inf = 10 −3 .The dotted red line and the solid blue line depict the numerical solutions of (20) and (19), respectively, with R given by (21) and keeping the same data.We can observe the transition from an initial stage in which increase in biofilm height dominates to a stage with faster horizontal spread.The green line is a reference self-similar approximation.

Figure 6 :
Figure 6: Layered distribution of dead, active, and inert cells, as illustrated by slices of a growing biofilm.Dead cells appear at the bottom of three peaks present in the initial biofilm seed.
producer with probability p e = cs K * s +cs 1 − cn K * n +cn .Cells become become inert with probability p i = 1 − cn K * n +cn .A non-surfactin-producer whose volume has surpassed a critical volume for division, divides with probability p d = cn K * n +cn .Figure 6 represents the cell type distribution for a growing biofilm.

Figure 7 :
Figure 7: Top view of the evolution of a biofilm seed formed by two layers of cells with a diameter of 40 cells, depicted at steps: (a) 45; (b) 90; (c) 135; (d) 180; (e) 225; (f) 270.Darkest colors correspond to layers of increasing height up to 10 cells.Contour undulations develop as the biofilm spreads.

.
The autoinducer concentrations are governed by balance equations of the form (16) with sourcesr surf = k surf (1 − c surf c surf +K surf ) and r cx = k cx (1 − ccx ccx+Kcx ), respectively, as well as no-flux boundary conditions at the biofilm boundaries.