## Abstract

In this contribution, we present a constitutive model to describe the mechanical behaviour of microbial biofilms based on classical approaches in the continuum theory of polymer networks. Although the model is particularly developed for the well-studied biofilms formed by mucoid *Pseudomonas aeruginosa* strains, it could easily be adapted to other biofilms. The basic assumption behind the model is that the network of extracellular polymeric substances can be described as a superposition of worm-like chain networks, each connected by transient junctions of a certain lifetime. Several models that were applied to biofilms previously are included in the presented approach as special cases, and for small shear strains, the governing equations are those of four parallel Maxwell elements. Rheological data given in the literature are very adequately captured by the proposed model, and the simulated response for a series of compression tests at large strains is in good qualitative agreement with reported experimental behaviour.

## 1. Introduction

Microbial populations in a matrix of extracellular polymeric substances (EPSs) adhering to surfaces and interfaces are generally referred to as biofilms [1–3]. Biofilms fulfil important ecological functions in nature and have various industrial applications. However, they are also responsible for many technical problems and associated with a number of diseases and infections [1,4]. The mechanical properties represent one of the many facets to be investigated for a better understanding of these complex systems, whose main constituents are microbial cells, water and EPS. The specific composition of the EPS depends not only on the micro-organisms but also on the particular environmental conditions under which it is produced and thus varies greatly [3]. Generally, however, the biopolymeric components comprise polysaccharides, proteins, lipids and nucleic acids [3,5–7]. Some of the EPSs are cross-linked into a network by permanent or temporary junctions. Upon hydration, this yields a hydrogel with only a few per cent of organic matter and a high water content, which is reported in a range of 87–99% [8–10].

Mathematical modelling concerns various fields in biofilm research, growth kinetics probably being the most active one of these; the interested reader is referred to previous reviews [11,12]. Focusing on the mechanical behaviour of biofilms, appropriate material models are indispensable to simulate and predict their response to external loads as they may appear in many biofilm-related technical processes and problems. This also includes modelling of biofilm growth, which may strongly be influenced by the mechanical environment, for example, spatial constraints or hydrodynamic conditions. Additionally, models are used to interpret and understand experimental results and thus complement the mechanical characterization procedures. The simplest modelling approaches for biofilms make use of empirical elastic or shear moduli to express a linear relationship between scalar stress and strain measures as in earlier studies [13–17]. The viscoelastic response has classically been interpreted in terms of linear spring–dashpot analogies such as the standard linear solid, the generalized Maxwell or Burgers model [18–24]. Expressing the relaxation function in terms of a Prony series, Towler *et al*. [25] used the latter model within a commercial finite element program to perform two-dimensional simulations. The nonlinear elastic properties at large strains were taken into account by neo-Hooke, Arruda–Boyce or St Venant–Kirchhoff hyperelastic models [26–28]. In several models for biofilm development, biofilm was treated as a viscous [29,30] or viscoelastic fluid [31]. The viscoelastic contribution was, however, often assumed negligible when time scales of growth much greater than the characteristic relaxation times were considered [32,33]. To include effects on a bacterial level, Alpkvist *et al*. [34] proposed a ‘hybrid’ model as combination of the individual-based modelling concept [35] and a viscous fluid description of the surrounding EPS matrix. Cogan & Keener [36] modelled biofilm as a gel consisting of EPS and water. Again, on the time scale of growth, the elastic stresses in the EPS were assumed to be relaxed, and the network was considered to be a Newtonian fluid, giving rise to a ‘two-fluid model’ [37]. An alternative ‘one-fluid two-component’ formulation was elaborated in the phase-field models by Zhang *et al*. [38,39], where the two components represent EPS and bacteria, on the one hand, and nutrients and solvent on the other. As a part of detailed studies of biofilm development, they investigated different constitutive equations for the EPS network, including two viscous and two viscoelastic models. Alpkvist & Klapper [40] simulated the mechanical response of biofilm, including detachment under fluid shear based on the immersed boundary method. The biofilm was discretized into a set of breakable springs, whose behaviour was assumed linear but could be extended to ‘more sophisticated’ behaviour.

In this paper, a viscoelastic model for biofilm is established and, particularly, biofilms produced by mucoid *Pseudomonas aeruginosa* strains are considered as a representative example. Although the structure and constitution of biofilm are generally complex, the EPS and particularly the exopolysaccharides seem to play the major role for mechanical stability and behaviour (cf. [3]). In line with other works in the literature, we consider the EPS matrix as a hydrated polymer network. However, instead of assuming *a priori* elastic or viscoelastic constitutive equations for this matrix, we allow for the properties of the polymer chains and their junctions. More precisely, the model network consists of worm-like polysaccharide chains linked by temporary junctions of different qualities. The polysaccharide alginate is the major carbohydrate component of the considered *P. aeruginosa* EPSs, which have largely been studied [3,10,41]. Alginates are produced by brown algae and bacteria, and represent linear copolymers of the monosaccharides *α*-l-guluronic acid (G) and *β*-d-mannuronic acid (M) [42]. Bacterial alginates differ from their algal counterparts by partial acetylation at the M residues [43,44] and, moreover, those produced by mucoid *P. aeruginosa* strains lack homopolymeric G-blocks [43,45,46]. Despite the absence of these blocks, which are usually associated with the formation of strong ion-mediated junction zones [47,48], the mechanical properties of *P. aeruginosa* biofilms are very sensitive to the concentration of divalent cations, particularly calcium (Ca^{2+}) [22,49]. This particular behaviour as well as the mechanical characteristics of *P. aeruginosa* biofilms and their EPS in general have been studied by various methods: Stoodley and co-workers [13,31,50,51] established a testing technique in which the deformation of a biofilm structure, a ‘streamer’, is measured while it is subjected to different fluid shear conditions. These studies reveal that the biofilm response is nonlinear and characterized by hysteresis in a loading and unloading cycle. Plotting wall shear stress acting on the streamer against engineering strain yielded J-shaped curves [13,31,51]. Körstgens *et al*. [14,18,49] accomplished compression experiments with a film rheometer and also observed strain stiffening. However, they attributed this nonlinear behaviour to the alignment of the upper platen of the compression device with the biofilm sample and not to the material properties. Initially, the biofilm stiffness increased with the strain rate, a decrease at higher rates was attributed to experimental limitations. When Ca^{2+} ions were added to the medium above a critical concentration, the biofilm became significantly stiffer. Moreover, biofilms without Ca^{2+} lost their integrity while swelling and revealed a quicker decay of stress in relaxation experiments [18]. The studies by Wloka *et al*. [21,22] focused on the rheological properties obtained in frequency and strain sweep tests with a rotational rheometer. They confirmed that the Cox–Merz rule holds, indicating the presence of entanglements in the hydrogel and determined changed rheological properties when calcium was added to the medium [22]. For a comprehensive review of available testing methods and mechanical characteristics of biofilms, the reader is referred to Böl *et al*. [52].

The remainder of this paper is organized as follows. In §2, we present the modelling of the biofilm as a viscoelastic polymer network. Numerical examples and a comparison with experimental data in the literature are given in §3. The results are discussed in §4 and, finally, conclusions are provided in §5.

## 2. A continuum description of the extracellular polymeric substance network

Taking the experimental findings into consideration, we establish a viscoelastic constitutive model for the EPS network of *P. aeruginosa* biofilms within the classical framework of nonlinear continuum mechanics and, particularly, methods developed for the constitutive modelling of polymer networks will be exploited. We will consider the alginate macromolecules as worm-like chains (WLCs) [53], and the eight-chain model [54] is used to represent a three-dimensional WLC network. At first, the network will be considered hyperelastic and after deriving the strain energy function, we address two important aspects of the model: its behaviour when the chain characteristics approach the Gaussian limit and its convexity properties. Finally, viscoelasticity is included by considering the network junctions to be transient.

### 2.1. A network of worm-like chains

The WLC is characterized by its contour length *L*, i.e. its fully extended length (figure 1), and the persistence length , which gives the typical length scale over which the directional correlation along the chain is lost [55]. Comparing these two lengths, one can distinguish flexible WLCs for which is much smaller than *L*. For these, an interpolation for the relation between the force *f* applied on the chain ends and its end-to-end distance *r* (figure 1) was suggested in Bustamante *et al*. [56] and was supplemented by a seventh-order polynomial correction by Bouchiat *et al*. [57]. Restricting this correction to a single square-term, Ogden *et al*. [58] finally obtained the interpolation formula
2.1where *k*_{B} = 1.381×10^{−23} J K^{−1} denotes Boltzmann's constant and *Θ* stands for the absolute temperature. An elegant way to relate chain and network behaviour was proposed by Arruda & Boyce [54], where eight chains spanning from the centre to the eight vertices of a cube are considered to represent the reference network. While the original eight-chain model was furnished with randomly jointed chains expressed in terms of Langevin statistics [59,60], modifications using Marko & Siggia [56,61] or MacKintosh [62] representations for WLCs were proposed in [63,64] and [65], respectively. The eight-chain representation leads to a simple expression for the non-affine stretch of the eight chains in terms of the first principal invariant of the right Cauchy–Green tensor *I*_{1} = tr **C** as [54]
2.2where *r*_{0} denotes the initial end-to-end distance in the reference state (figure 1). By introducing the ratios *Λ* = *r*_{0}/*L* and , integrating the force (2.1) and inserting relation (2.2) yields the strain energy function of the corresponding eight-chain model as
2.3where the factor *n* represents the number of chains per unit reference volume and *W*_{0} is a constant term that guarantees zero strain energy in the reference state. Note that here and throughout this section, the chains will be assumed to be inextensible and the network incompressible, which implies for the determinant of the deformation gradient **F** that det **F** = 1 and 3 ≤ *I*_{1} < 3/*Λ*^{2}.

For very long chains , one would expect that, in the limit, the network acts like a Gaussian one and that, consequently, the strain energy recovers the neo-Hookean form. Ogden *et al*. [66] have recently highlighted that eight-WLC models of the form (2.3) fail to reproduce this state. This problem is remedied by considering that the mean squared end-to-end distance *r*_{0}^{2} is related to the contour and persistence lengths by [55]
2.4

If , one finds the limit and thus and . Equation (2.3) then takes the form 2.5

Bearing in mind that , we consider the series expansion of the first term in the curly brackets around . Because , the series converges [67], and we write 2.6

Inserting this result into equation (2.5), we obtain the strain energy function 2.7

Accordingly, the strain energy function (2.7) differs from the neo-Hookean model by the constant term in brackets, and the remainder of the series expression. Provided that the latter is negligible for moderate values of . Consequently, the deviation from the neo-Hookean response will be very small for relevant states of deformation. This is illustrated in figure 2, where the strain energy according to equation (2.3) is plotted against for different ratios of in units of in comparison with the neo-Hookean model. Apart from being physically reasonable, this result also provides the link between the current approach and several models that were applied to biofilms before, as will be discussed later.

For the applicability of hyperelastic models, the convexity properties play a critical role. For a class of isotropic and transversely isotropic eight-WLC models, these properties were studied by Kuhl *et al*. [68] in uniaxial and biaxial tension as well as pure and simple shear. We note that the strain energy function given in (2.3) fits into the framework of polyconvexity [69], which applies to generic load cases. A proof is provided in appendix A.

### 2.2. Presence of different junction types

A chain is usually considered a sequence of monomeric units connecting two cross-links [60] and can thus be substantially shorter than the macromolecule itself. The cross-links are created by physico-chemical interactions between the macromolecules and may be of different type [3,6]. Potential interactions responsible for cross-linking the EPS matrix are electrostatic interactions such as those mediated by divalent ions, dispersion forces and hydrogen bonds [10]. Moreover, as a fourth possible mechanism, theoretical considerations predict and rheological data indicate the presence of entanglements in the EPS matrix from *P. aeruginosa* [14,21,22,49]. On the contrary, covalent bonds seem to play no significant role [6,70], in clear contrast to rubber networks. We assume that all four former connection mechanisms are present with different amounts, suggesting an additive decomposition of the chain density
2.8into four parts *n*_{i}, *i* = 1,2,3,4, representing the number of chains per reference volume connected by dispersion forces, entanglements, hydrogen bonds (*i* = 1,2,3)^{1} and electrostatic interactions (*i* = 4), respectively. Intuitively, with the number of junctions, the chain length *L*_{i} and thus the mean distance *r*_{0i} between two junctions of type *i* will vary according to equation (2.4). We derive the strain energy of a network with four different types of junctions from equation (2.3) as
2.9where is again a constant term that guarantees zero strain energy in the reference state and
2.10with and . According to classical principles in continuum mechanics, the second Piola–Kirchhoff stress tensor **S** results in consideration of the incompressibility constraint as [71]
2.11where **I** is the identity tensor of second order, *p* denotes the arbitrary hydrostatic pressure and represent response functions in terms of the first principal invariant.

### 2.3. Temporary junctions

In the preceding sections, the cross-links were assumed to be permanent, and equation (2.11) yields the stress in an elastic network. The experimentally reported viscoelastic behaviour, however, can be interpreted in terms of a temporary network in which the junctions are ‘fluctuating’ (cf. [14,49]). Consequently, relaxation and creep are enabled by perpetual releasing and reconnecting of chains. These considerations are in line with Green and Tobolsky's theory of relaxing media [72], according to which the rate of breaking of any set of chains is given by the simple differential equation
2.12with *T*_{i} being a constant that characterizes the ‘lifetime’ of a junction. According to the theory [72], the overall number of chains remains constant and whenever one chain breaks, another one is reformed in an unstrained state. The deformation experienced by a chain depends thus not only on the current time *t* but also on the point in time *s* ≤ *t* at which it was formed. This deformation is expressed by the relative deformation gradient between the times *s* and *t* [72]
2.13

With this and the definition at hand, one finds the invariant
2.14where the colon indicates the scalar product of two second-order tensors. By analogy with equation (2.2), represents the current relative non-affine stretch of chains that have formed at time *s*. Because the initial number *n*_{i} of such chains decreases with exp(−(*t* − *s*)/*T*_{i}) according to the solution of equation (2.12), their contribution to the stress at time *t* calculates as
2.15

Because the number of chains reformed during d*s* and surviving until *t* is (cf. [72]), the stress **S**_{i} at time *t* results from integration over the complete history as
2.16where the memory function has been introduced. The integral expression (2.16) represents a generalization of Green and Tobolsky's theory to nonlinear material behaviour and takes the form of a particular K-BKZ equation [73–75] in Lagrangian description. Finally, from equation (2.11), we deduce the Cauchy stress tensor for the viscoelastic material as^{2}
2.17where we bear in mind that depends on the relative deformation of the network between times *t* and *s*.

## 3. Comparison with experiments and numerical examples

In this section, the viscoelastic model is compared with experimental results. At first, the linearized case, valid for small strains, is considered and applied to estimate material parameters from rheological data. On the basis of these parameters and additional assumptions, large strain compression experiments are simulated, which are motivated by different compressive testing procedures in the literature.

### 3.1. Small strain shear response

As a consequence of the behaviour illustrated in figure 2, Green and Tobolsky's model is approached for up to a scalar factor (cf. [72, p. 90]). Similar implications hold for small strains and we particularly consider the shear stress response. In simple shear, the deformation gradient takes the form , where *γ*(*s*) denotes the time-dependent amount of shear. Restricting to small strains, , one finds and neglecting terms of higher order, the stress (2.17) reduces to
3.1where results from the response functions evaluated at the limit . Accordingly, the shear stress for small strains reads
3.2

The corresponding rate equation is obtained by taking the time derivative on both sides
3.3and recovers the form of a rheological model that consists of four parallel Maxwell elements. The frequency-dependent storage and loss moduli characterizing the response to small oscillatory shear are thus given by [75,76]
3.4where *ω* denotes radian frequency. These relations are particularly useful to estimate the time constants *T*_{i} from standard rheological experiments. To this end, we consider the results by Wloka *et al*. [22], who tested biofilms formed by several strains of *P. aeruginosa* with a parallel plate rheometer. Their study included biofilms of the strain FRD1 grown with either 1 or 10 mmol l^{−1} CaCl_{2} added to the medium. In order to estimate the values of *T*_{i} and *μ*_{i}, the least-squares error between the frequency-dependent storage and loss moduli provided in Wloka *et al*. [22] and the expressions (3.4) was minimized with Maple software (Maple 15, Maplesoft). According to our assumptions, the times *T*_{i} should be characteristic of the type of junction and should thus not vary with the amount of Ca^{2+}, in contrast to *μ*_{i}, which relates to the number of junctions *n*_{i}, the contour lengths *L*_{i} and the persistence length . Thus, we determined four time constants and two sets *μ*_{i}^{1} and *μ*_{i}^{10} referring to 1 or 10 mmol l^{−1} CaCl_{2}, respectively. The results and corresponding parameters are given in figure 3 and table 1, respectively. Note that the ratios of the moduli *μ*_{i}^{10}/*μ*_{i}^{1} reveal a more than 150-fold increase for the Maxwell-branch with the largest relaxation time (*i* = 4), which, in the model, is associated with chains connected by Ca^{2+} junctions. The ratio for the other branches increases as well, even though by a factor less than 5.

### 3.2. Large strain compressive response

Compression experiments represent a relatively straightforward method for characterization of biofilm behaviour at large strains. Compression of the biofilm has been accomplished between either two parallel platens [14,18,49] or by pushing an indenter into the biofilm [77,78]. Although the samples are expected to be highly constrained and deformations inhomogeneous in the latter cases, we illustrate the constitutive equations for homogeneous uniaxial compression with lateral extension because this simple mode serves well to highlight the basic characteristics of the model. Given the stretch *λ*(*t*) as the ratio of current to initial height of the biofilm at time *t*, the deformation gradient takes the form
3.5

Assuming that the lateral directions are free of traction, the stress difference *σ*_{11} − *σ*_{22} serves to eliminate the hydrostatic pressure. From equation (2.17), we thus obtain the compressive component of the Cauchy and first Piola–Kirchhoff stress, respectively, as
3.6

In the following, we will illustrate the compressive stress response for several loading conditions, where Maple software was used to evaluate the integral in equation (2.17) numerically.

#### 3.2.1. Material parameter estimation

While the time constants *T*_{i} can directly be identified from small strain rheological data, there is a non-unique dependence of the parameters *μ*_{i} on *L*_{i}, *n*_{i} and for each *i* and both Ca^{2+} concentrations. The results in table 1 provide only eight equations for the 18 unknowns , , as
3.7where the superscripts again indicate the values associated with the different amounts of Ca^{2+} in the medium. Generally, calculating the least-squares error between model predictions and suitable experimental data that reveal the nonlinear behaviour of the material at large strains could be used to estimate these parameters. However, in lack of such data, we use the following assumptions, which will be discussed in §4: assuming that is independent of the Ca^{2+} concentration, we set nm (cf. [79]). Moreover, Wloka *et al*. [22] documented the actual calcium concentrations *c*_{1} and *c*_{10} in the biofilm, which were 0.047 and 0.334 mmol l^{−1} and thus significantly less than the concentrations in the medium. Having the number of Ca^{2+} ions per volume, we assume that five ions are necessary to form a stable junction and that these junctions are distributed uniformly in space. With this, a rough estimate of the mean distance between two junctions can be obtained from the distance of nearest neighbours (cf. [80]) as , *β* = 1,10, where *N*_{A} = 6.022 × 10^{23} mol^{−1} is Avogadro's constant. Inserting this into equation (2.4), we obtain the estimate
3.8which, for a given , can be solved for corresponding to the two Ca^{2+} concentrations. The values for *n*_{4} then follow implicitly from equation (3.7). The remaining contour lengths are set to *L*_{i}^{1} = *L*_{i}^{10} = 200 nm, *i* = 1,2,3, providing the missing six equations, and all calculations were referred to *Θ* = 294 K. The parameters defined in this way and used throughout this section are given in table 2.

#### 3.2.2. Rate dependence

At first, we assume that the specimen is compressed with a constant rate *α*_{λ}, say, so that *λ*(*t*) = 1 − *α*_{λ}*t*. Varying *α*_{λ} over four orders of magnitude from 0.0001 s^{−1} to 0.1 s^{−1}, the responses of both the low and high Ca^{2+}-biofilm increase with rate (figure 4), where the stress–stretch curves at rates of 0.0001 s^{−1} and 0.001 s^{−1} are nearly indistinguishable for the former.

#### 3.2.3. Relaxation

Next, we simulate relaxation experiments, which have widely been used to characterize the viscoelastic behaviour of biofilms [52]. In the compressive mode, the biofilm is thereby rapidly compressed to a certain level and then kept at this thickness while monitoring the stress response. In this case, the stretch is given by 3.9

Here, we use *α*_{λ} = 1 s^{−1} and simulate biofilm compression to either *λ* = 0.95, *λ* = 0.75 or *λ* = 0.5. The normalized stress *P*_{c}(*t*)/*P*_{c}(*t*_{1}) against time is plotted in figure 5 and reveals the distinct differences in the relaxation behaviour of the two biofilms over a period of 200 s. The underlying relaxation mechanism is once again depicted in figure 5, which provides illustrations of the network shortly after loading and in a more relaxed state: owing to the perpetual chain breakage and reformation, the load-bearing chains diminish and are replaced by new, unstretched ones. The longer relaxation times experimentally observed for biofilms grown under higher Ca^{2+} conditions [18] are thus explained by the greater amount of chains with the long relaxation time *T*_{4}.

We note that the relaxation experiment provides a straightforward analytic solution for the stress if the ramp loading in equation (3.9) is idealized as a step from *λ* = 1 to *λ* = *λ*_{1} at *t*_{1}. In this case, the compressive stress is obtained as (see appendix A)
3.10where , and is characterized by linear viscoelastic relaxation behaviour. For *λ*_{1} = 0.5, this is illustrated in figure 6 as a double logarithmic plot over 5000 s. From this representation, the four relaxation modes associated with the four relaxation times *T*_{i} can easily be identified and are indicated by dashed lines.

#### 3.2.4. Influence of swelling

The water content of the hydrated EPS network depends on the osmotic pressure, which is balanced by the elastic response of the stretched network (see [81] and references therein). Because alginates are polyelectrolytes [42], the osmotic pressure depends not only on solvent quality, i.e. the free energy of dilution, but also contains an ionic contribution driven by the differences of the chemical potential inside and outside the network (cf. [81,82]). In the scope of the present contribution, we omit a detailed discussion of the swelling behaviour but take into account that changing the ion concentration in the surrounding medium can cause a volumetric change of the network. Denoting the swelling volume ratio by and assuming isotropic behaviour, the volumetric deformations can be expressed as . The deformation gradient **F** with respect to the network in the state as it was formed and the according first principal invariant are then given by
3.11where represents the deformations applied to the swollen network (cf. [83]). Note that, as a consequence of relation (2.2), may also be understood as the change of initial end-to-end distance of two junctions, which is equal to in the swollen state as illustrated by the inset in figure 7.

For a transient network, swelling competes with the relaxation processes, which may prevent an equilibrium between osmotic pressure and network stress. In order to shed light on the special behaviour of swollen biofilms observed by Körstgens [18] and discussed in §1, we assume that the swelling is defined by a constant rate such that after a time of the volume doubles. We simulate a biofilm, which starts swelling at *t* = 0 and from *t* = *t*_{1}, it additionally undergoes compression with a constant rate of *α*_{λ} = 0.005 s^{−1}. This is expressed by the deformation gradient
3.12with
3.13where we chose *t*_{1} = 10 min and *α*_{ν} = 0.1 min^{−1}. Note that the expressions for the stress (2.17) and (3.6)_{1} are still valid; however, the constraint of mechanical incompressibility obviously takes the form in this case. We refer the measured force to the geometric state at *t* = *t*_{1} when the compression test started. The nominal stress tensor relating to this configuration calculates as
3.14and for the compressive component, one obtains
3.15

The strikingly different behaviour of the two biofilms upon swelling is highlighted in figure 7. While the Ca^{2+}-rich biofilm stiffens significantly when compressed in the swollen state, the low Ca^{2+} biofilm becomes weaker. This behaviour is in qualitative agreement with the experimental results provided in earlier studies [18,70] for *P. aeruginosa* strains SG81 and FRD1. In the model, the successive expansion of the network by the entering water is relaxed by the transient nature of the junctions. Thus, when the rate of expansion is higher than the rate of chain breakage and formation, the network response becomes stiffer. In the given example, this is the case for the Ca^{2+}-mediated junctions, which explains the stiffer response.

## 4. Discussion

The model proposed in §2 provides a viscoelastic extension of the Arruda–Boyce eight-chain network model furnished with WLCs. The additive decomposition of the network contributions in equation (2.9) states that the response of one network with four different types of junctions is equal to the summed response of four separate networks with one type of junction, respectively. This simplified treatment of the network provides some important advantages for the model formulation. For example, the limiting behaviour and convexity properties discussed in §2.1 are maintained by superposition and, finally, the characteristics of four parallel Maxwell elements are captured for small shear strains. The comparison of the storage and loss shear moduli resulting from the linearized model provided in figure 3 with rheological data on a *P. aeruginosa* FRD1 biofilm [22] shows excellent agreement. This indicates that the four discrete relaxation times that result from the network decomposition are sufficient to characterize the rheological behaviour of the biofilm within the experimentally studied frequency domain between 0.01 and 100 rad s^{−1}. Notwithstanding, additional terms may be included to capture possible effects beyond these limits. It is worth mentioning that the largest identified relaxation time of 1925 s 32 min is of the same order of magnitude but around 1.8-times larger than the value of 18 min obtained by Shaw *et al.* [84], as a common value for a wide range of biofilm samples.

The techniques for characterization of mechanical biofilm properties at large strains are little standardized (cf. [52]), the exact boundary conditions often unclear and the deformations highly inhomogeneous. This renders the quantitative interpretation of experimental results and a comparison with constitutive models difficult. For this reason, the estimation of material parameters for the finite strain model was based on a few assumptions. While the contour lengths of alginate chains between junctions of type 1–3 were set to 200 nm for illustration purposes, the remaining values were motivated by experimental findings. Particularly, we set the persistence length constant to 15 nm according to a study by Vold *et al*. [79], who determined for several alginates of different origin at a fixed ionic strength and we note that values in a similar range of 11.5 nm for M-rich to 26 nm for G-rich alginates were discussed by Gamini *et al*. [85], and references therein. However, owing to the polyelectrolyte character of the alginates, may depend on the ion concentration (cf. [42]) and thus change with an increased amount of calcium in the biofilm. The average distance between nearest Ca^{2+} junctions was estimated from measured ion concentrations in the biofilm and identified with the initial end-to-end distance of a chain. Thereby, we assumed that a stable junction is formed by five ions motivated by the junction zones of the egg-box model (cf. [47,48]). Although alginates produced by mucoid *P. aeruginosa* strains lack homopolymeric G-blocks [43,45,46], which are a prerequisite for the classical egg-box model, junction zone formation involving alternating sequences has been suggested [86], which may play a role for the gelation in bacterial alginates [87].

Even though based on assumptions, the calculated parameters have some interesting implications. Just as the ratios, the parameters for the high Ca^{2+} biofilm differ not only for the Ca^{2+} network but also for *i* = 1,2,3. The sum is a measure of the total polymer content in the network. The parameters given in table 2 yield indicating a strong increase in alginate per volume. This is in agreement with a study by Sarkisova *et al*. [88], who reported that an increase in Ca^{2+} leads to strong changes in structure and composition of the extracellular matrix of FRD1 biofilms as well as an increase in alginate content. Moreover, the overall number of chains given in table 2, , is in agreement with that of ‘elastically effective chains’ estimated as in Wloka *et al*. [22].

The simulation results provided in §3.2 are generally in good agreement with experimental results such as those given in Körstgens [18]. The model predicts stiffening with increasing strain rate (figure 4). The high Ca^{2+} biofilm shows a throughout stiffer response and stress relaxation takes significantly longer (figure 5). Finally, the simplified treatment of time-dependent swelling revealed the ability of the model to capture the distinct difference in mechanical behaviour between high and low Ca^{2+} biofilms subjected to swelling (figure 7).

Putting the suggested model into perspective with previous approaches reveals that many of these are contained as special cases. At the end of §2.2, we have shown that for long chains, the neo-Hookean rubber model is recovered from the proposed eight-chain approach. This model was used by Dupin *et al*. [26] for the simulation of biofilms and also by Zhang *et al*. [38] as one limit case in a larger study. The latter authors also applied the Lodge rubber-like liquid or upper convected Maxwell model (cf. [75,76]), which is the Gaussian limit case of the transient network model in equation (2.17). The original Arruda–Boyce strain energy function [54] was used by Böl *et al*. [27] and is conceptionally similar to the proposed model in the elastic limit given in equation (2.3). Finally, the generalized Maxwell model with material functions (3.4) is a standard approach to describe the viscoelastic behaviour of *P. aeruginosa* biofilms in the frequency [21,22] and time domain [18]. In this sense, the model presents a valuable advancement because it provides a fully three-dimensional and nonlinear continuum mechanical description of the EPS network while still being able to recover many previous approaches.

Nevertheless, there are some limitations to be addressed in future work. First of all, the biofilm properties have entirely been assigned to the alginate network. This main responsibility of the polymer network and, potentially, the solvent for biofilm mechanics was assumed by others before (cf. [26,36,38,39]) but it was also noted that other EPS constituents and cells may contribute [22, p. 282]. We also remark critically that with the simple transient network approach presented herein, the behaviour predicted by the Cox–Merz rule cannot be captured in contrast to what is seen in experiments [22]. Consideration of more elaborate evolution equations for the number of chains may serve to remedy this deficiency. Finally, the framework presented herein focuses on the elastic and viscoelastic behaviour. In its current form, it neglects plastic yielding, which was considered to interpret some experimental results [14,49].

In order to study and highlight some basic features of the proposed model, the simple case of homogeneous uniaxial compression with lateral extension was considered. A numerical treatment of the constitutive equations including the implementation in a finite element framework and more demanding computations will be provided elsewhere.

## 5. Conclusions

In the present contribution, we have proposed a constitutive theory for microbial biofilms based on the premise that the EPS is mainly responsible for the mechanical behaviour and can be represented by a network of polysaccharide chains. The model was particularly developed for biofilms formed by mucoid *P. aeruginosa* strains as a result of the large interest they have gained in the literature. The framework could, however, easily be adapted to other biofilms as it captures many of the characteristics generally observed in biofilm mechanics. The model consists of four superposed WLC networks, each connected by transient junctions of a certain lifetime. The resulting constitutive equations take the form of a particular K-BKZ type fluid and include various previous models applied to biofilms as special cases. Particularly, for small shear strains a four-branch generalized Maxwell model is recovered, which favourably compares with experimental rheological data. Fitting the model to these data immediately yields the time constants and additionally a set of implicit equations to be satisfied by the remaining material parameters, which are the number of chains in each network, their contour lengths and the persistence length of the alginate. For illustration purposes, further assumptions were made to estimate these parameters and the model was applied to simulate compression tests. The results are in general agreement with experiments and can particularly explain some of the striking mechanical differences between biofilms grown at different concentrations of calcium.

## Acknowledgements

Partial support for this research was provided by the Deutsche Forschungsgemeinschaft under grant no. BO 3091/6-1. The authors are grateful to A. Bolea Albero for valuable discussions.

## Appendix A

**A.1. Polyconvexity of the strain energy function**

Polyconvexity of the strain energy function *W* implies the existence of a convex function *W* such that [69,89], being the cofactor of **F**. The function (2.3) depends exclusively on the first principal invariant and thus we prove convexity with respect to **F**. Because *I*_{1} is a convex function of **F** (cf. [90]), the strain energy function (2.3) is polyconvex if is convex and monotonically increasing according to the monotone composition of convex functions (cf. [90, lemma B.9]). Thus, it is sufficient to ensure the first and second derivatives of *W* being positive. After some algebra, we calculate these derivatives as
A1and
A2

Recalling that 0 ≤ *Λ* < 1 and *I*_{1} < 3/*Λ*^{2}, the first two terms in expression (A 2) are non-negative. Moreover, the third term is easily shown to have a minimum of 1.25 at *I*_{1} = 6.75/*Λ*^{2}. Thus,
A3implies . Consequently, *W* is a convex and moreover a monotonically increasing function of *I*_{1}. Accordingly, we infer that
A4which implies monotonicity of *W*.

**A.2. Analytic solution of stress relaxation**

For step loading, for and for *t* ≥ *t*_{1}. Thus, in both intervals, **F** is constant and with *J* = 1, the integral (2.17) can be separated as
where is the left Cauchy–Green tensor. Inserting equation (3.5) at time *t*_{1}, we obtain for the first stress difference
A5Dividing by *λ*_{1} yields the nominal stress as given in equation (3.10).

## Footnotes

- Received August 22, 2012.
- Accepted September 10, 2012.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.