Network patterns in exponentially growing 2D biofilms
Cameron Zachreson, Xinhui Yap, Erin S. Gloag, Raz Shimoni, Cynthia B. Whitchurch, Milos Toth
NNetwork patterns in exponentially growing 2D biofilms
Cameron Zachreson, Xinhui Yap, Erin S. Gloag,
2, 3
Raz Shimoni, Cynthia B. Whitchurch, and Milos Toth School of Mathematical and Physical Sciences, University of Technology Sydney, Ultimo, NSW, 2007, Australia The ithree institute, University of Technology Sydney, Ultimo, NSW, 2007, Australia Center for Microbial Interface Biology, Ohio State Univeristy, 460 W. 12th Ave., Columbus, OH 43210 (Dated: November 6, 2018)Anisotropic collective patterns occur frequently in the morphogenesis of 2D biofilms. These pat-terns are often attributed to growth regulation mechanisms and differentiation based on gradients ofdiffusing nutrients and signalling molecules. Here, we employ a model of bacterial growth dynamicsto show that even in the absence of growth regulation or differentiation, confinement by an enclosingmedium such as agar can itself lead to stable pattern formation over time scales that are employedin experiments. The underlying mechanism relies on path formation through physical deformationof the enclosing environment.
PACS numbers: 87.18.Fx, 87.17.Jj, 87.18.Gh
In surface-associated bacterial colonies, growth regu-lation caused by a diffusion-limited nutrient supply hasbeen suggested as a primary mechanism of pattern for-mation, resulting in branched, locally anisotropic mor-phologies [1, 2]. However, diffusion-limited growth doesnot tell the whole story. Mechanical processes are knownto play essential roles in morphogenesis [3], and we showhere that deformation of the enclosing environment bymoving bacteria can lead to the emergence of persis-tent, stable, anisotropic network patterns such as those inFig. 1(a,b) even under conditions of exponential colonygrowth. Furthermore, we present strong evidence thatin
Pseudomonas aeruginosa , well-studied trail follow-ing mechanisms involving the deposition of extracellularpolymeric substances (EPSs, or “slime”) cannot producethese patterns in the absence of such mechanical effects.Interstitial biofilms form when bacteria are introducedto the space between two apposed surfaces and are con-strained to a quasi-2D arrangement. If one surface is ahydrogel such as agar or gellan gum [4, 6], the bacteriahave access to nutrients diffusing through this medium[7]. Here we consider such a situation, where nutrientdepletion is negligible and the bacterial growth rate isnot retarded in the colony interior. This assertion is sup-ported by histograms of individual growth rates we ob-tained from of the edge and interior of an expanding
P.aerugionsa interstitial biofilm [Fig. 1(a,b), Supplemen-tal Movies S1, S2 [5] respectively]. The histograms showthat growth rates measured in the colony interior and thecolony edge are approximately equivalent [Fig. 1(c)], andnutrient depletion therefore cannot be used to explain thenetwork morphology seen in Fig. 1(a,b), characterized bydense trails of bacteria and voids of zero density.Previously, Gloag et al. observed that these bacte-ria can deform the soft enclosing material through whichthey move, and confinement in the resulting ‘furrows’appears to contribute to network pattern formation [4].However,
P. aeruginosa and many other surface-motilebacteria [8] move through the extension and retraction b) a) c) e) P( k g ) -3 k g edge interior model d) f) g) coverslip substratum interioredge FIG. 1. (a, b) Phase-contrast images of the edge and inte-rior (respectively) of an expanding
P. aeruginosa interstitialbiofilm [4]. Each image is frame 500 from a total of 1000 usedto make Supplemental Movies S1 and S2 [5]. Binary masks ofthe cells produced by the image segmentation algorithm areshown in red. (c) Histograms of the probability P for a cellto be growing at rate k g , as calculated from the movies cor-responding to the colony edge (a, solid bars in c) and interior(b, dashed bars in c). The dashed line represents the growthrate used in our simulations. (d-g) Schematic illustrations ofthe processes comprising our model of bacterial behavior. of type IV pili (T4P) [9], which have increased bindingaffinity to extracellular polymeric substances (EPS) se-creted by the moving bacteria [10], a phenomenon thathas also been suggested to result in trail-following be-havior [11, 12]. To investigate the relative importance ofthese two processes in the emergence of the observed net-work morphology, we carried out simulations of motilebacteria that interact with the hydrogel environment. a r X i v : . [ phy s i c s . b i o - ph ] J un d)e)f)a)b)c) FIG. 2. Snapshots of three biofilms simulated by the modelwith substratum stiffness coefficients of γ = 0 . , . , . . , . .
2% (d, e, and f respec-tively). Scale bars represent 100 w (a-c), or 100 µ m (d-f). The behavioral rules implemented in our simulations areillustrated in Fig. 1(d-g). Each simulation was initiatedwith a single motile, growing cell, and terminated after12 cell division cycles ( t f = 1 . × s ). Our resultsverify the importance of hydrogel properties in intersti-tial biofilm morphogenesis. By increasing the substratumstiffness in our model, colony morphology [Fig. 2(a-c)]was altered in qualitative agreement with experimentalresults obtained by increasing hydrogel monomer concen-tration [Fig. 2(d-f)]. [See [13] for details on the intersti-tial biofilm culture techniques used to produce the datain Figs. 1(a,b) and 2(d-f).]Bacterial behavior emerges from a complex interplayof many non-equilibrium and stochastic processes. Iso-lation of factors that are fundamental to morphogenesismechanisms requires a model that describes the complex-ity of the real system, with an experimentally constrainedparameter space. Here, we focus on the processes of cellgrowth and motility as the fundamental model ingredi-ents. The essential properties of the system that can-not be estimated a-priori relate to bacterial movementdynamics and individual growth rates. We therefore de-signed our biophysical model so as to enable incorpora-tion of experimentally determined time-scales into theseprocesses.The model used here is similar to one we reported pre-viously [14], except that it includes cell growth and di-vision. We parameterized our model using experimental data that quantify single-cell dynamics [9], and our ownmeasurements of individual growth rates [Fig. 1(c)] [5].Individual bacteria are simulated as self-propelled,capped rods of unit width w , that undergo repulsive colli-sions with one another [Fig. 1(d)], elongate at a constantrate k g and divide [Fig. 1(e)] at a critical length l max =7 w into two cells each with length l min = 3 w ± δl ,where δl represents asymmetric division and is randomlyselected at the time of division from the uniform interval δl ∈ [0 , w . Cell-cell collisions produce force and torqueaway from the contact point. It is also possible for cellsto pull on one another with their pili because we modelthe cell motility mechanism explicitly [14] by simulat-ing the extension, binding, and retraction of T4P [15].The equations for translational and rotational motion ofrods and our repulsive interaction scheme are the sameas those used in [2, 14, 16]. During T4P-mediated surfacemotility, P. aeruginosa cells spontaneously reverse polar-ity, this process is believed to be essential to chemotaxisphenomena, but polarity reversals occur stochasticallyeven in the absence of attractant/repellant gradients [17].Here, stochastic polarity reversals were implemented us-ing a ‘countdown’ type algorithm where the time betweenreversal events is selected from a Gaussian distributionwith a mean (cid:104) T rev (cid:105) = 1000 s and a standard deviation σ rev = 200 s .Substratum deformation [Fig. 1(f)] is simulated by asurface potential U that resists motion up its gradient. U increases when a bacterium is present, representinga depression in the substratum. Three parameters con-trol the mechanical properties of the substratum: thestiffness γ , and the deformation and restitution rates k U and β U of the potential U , respectively. The stiff-ness, γ = dU/dC s modulates the force (cid:126)F s exerted bysubstratum deformation, where C s represents surface to-pography, such that (cid:126)F s = − γ ∇ C s is the force exerteddue to deformation. The deformation rate per unit area, k s = k U /γ and restitution rate β s = β U /γ , define thesubstratum dynamics for a spatial region of area [∆ x ] so that ˙ C s = k s [1 − C s ( t ) /C max ][∆ x ] − β s C s ( t ), where C max = w and the accumulation term only applies if abacterium is present. These rate constants k s and β s areinversely proportional to γ and we hold their ratio fixedat β s /k s = 0 . U influences themotion of the bacteria.We simulate the function of secreted EPS in motility[Fig. 1(g)] with a trace left behind by bacteria that mod-ulates the T4P binding probability: P s = KP max + (1 − K ) P min , where K = C p / [∆ x ] is the fraction of localarea covered by binding sites for pili and C p is the localtrace value that accumulates and degrades analogouslyto C s according to the EPS deposition and degradationrates ( k eps = 0 . β eps = 5 × − , respectively).A similar model of motility bias due to EPS secre-tion was published recently [11, 12]. However, our modelexplicitly simulates the stochastic process of pilus bind-ing and retraction [14], instead of making a mean-fieldapproximation based on the assumption of a large num-ber of T4P per cell. Based on available literature, suchan assumption may not be valid for P. aeruginosa [9].This apparently subtle difference is significant becauseour findings suggest that the stochastic nature of T4Pbinding and retraction yields ineffective following of EPStrails that is not robust in the absence of furrowing. Thisis exemplified by Fig. 4(a), which illustrates a typicalcolony morphology in a case where physical interactionswith the substratum are weak. In such a case trail net-works could, in principle, emerge exclusively via the EPSstigmery mechanism but do not, making the EPS phe-nomenon alone unlikely to account for network forma-tion in
P. aeruginosa . A detailed discussion of this im-portant result, including tests of model variants lackingsubstratum deformation, is provided in the SupplementalInformation [5].Because the furrowing process is governed by the sub-stratum deformation rate k U and stiffness γ , these pa-rameters control the path formation rate and degree towhich cells are confined to the paths, respectively, mak-ing them critical to this trail formation mechanism. Toinvestigate their roles, we simulated bacterial growth un-der various environmental conditions defined by k U and γ , and recorded the colony configuration (cell positions,lengths, and orientations) at 20 s intervals.Three different morphological classes emerge as a func-tion of the parameters defining the mechanical proper-ties of the substratum. These consist of two qualita-tively distinct isotropic states, and a locally anisotropicnetwork similar to the one observed in experiments.Dense configurations become isotropic due to buckling[18] (Fig.3, γ = 1 . k U = 0 . γ = 0 . , . k U = 0 . γ = γ = γ = γ = k U = k U = FIG. 3. Plots of simulated furrow depth after 11 divisioncycles. Four different substratum stiffness coefficients ( γ ),and two different deformation rates ( k U ) are represented. Thescale bar represents 400 w . Unregulated exponential growth eventually outpacesall other processes in the model and acts to homoge-nize the colony morphology as t → ∞ . However, on thetime scales probed in typical experiments, the colony canrapidly become isotropic, [Fig. 4(a)] or it can grow in anetwork [Fig. 4(b)], depending on the values of k U and γ .The changes in colony morphology during growth seen inFig. 4(a,b) correlate with changes in the radial distribu-tion of cells around the colony center of mass [Fig. 4(c,d)],which transitions from a multi-modal to a uniform dis-tribution as the colony grows from a single cell movingback and forth in the hydrogel to become isotropic.The entropy of such radial distributions as a func-tion of length scale r and distance from the colony edge d edge provides a metric of local anisotropy in the densitydistribution, a signature of the network morphology ob-served in experiment [Fig. 1(a,b), and Fig. 2(d,e)]. Ouranisotropy parameter ∆ S φ ( r, d edge ) reflects the entropyof local distributions relative to that of the uniform dis-tribution. ∆ S φ takes a value of 1 for a 1D density distri-bution and approaches 0 for an isotropic (circular) dis-tribution [5]. P ( φ ) φ replication cycle -300-200-1000100200300 3002001000-100-200 k U = 0.001 γ = 0.5replication cycle
12 10 8 6 4 2 b) d) P ( φ ) φ replication cycle a) c) -300-200-1000100200300 3002001000-100-200-300 k U = 0.05 γ = 0.5replication cycle
12 10 8 6 4 2
FIG. 4. (a, b) Colony outline (defined as the area over which C s (cid:54) = 0) plotted after a number of cell replication cycles onsoft (a) and hard (b) agar (Supplemental Movies S6 and S4,respectively [5]). (c, d) Histograms of polar coordinates φ ofthe constituent bacteria integrated over the indicated replica-tion cycles for soft (c) and hard (d) agar. Quantitative comparison of ∆ S φ in the simulated andexperimental morphologies depicted in Fig. 2 are pro-vided in Fig. 5. The results indicate good agreement incases where network morphology is observed [Fig. 5(g,h)], with high degrees of anisotropy at all length scalesnear the colony edge, falling off to shorter length scalestowards the interior. For the dense morphology, quanti-tative similarity is confined to shorter length scales, closeto the colony edge (Fig. 5i). This discrepancy is likely dueto greater cell numbers and colony area in experiment,as well as an overestimation of the substratum resistancevalue γ in the corresponding simulation. Qualitatively,1 r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge adg b ce fh i r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge r d edge Δ S φ Δ S φ Δ S φ FIG. 5. Intensity maps of ∆ S φ as a function of length scale r and distance from the colony edge d edge . Measurements fromimages of bacterial biofilms grown on 0.5%, 0.8%, and 1.2%gellan gum substrata (a, b, and c respectively), are comparedto data measured from simulations with increasing γ valuesof 0.25, 0.5, and 1.0, (d, e, and f). Residuals comparing sim-ulation results to experimental measurements are shown in g,h, and i ( g = | a − d | , h = | b − e | , and i = | c − f | ). Alldistances are in units of average cell width w . however, there is an inversion in the scaling of ∆ S φ as afunction of r between the network and dense morpholo-gies in both simulation and experiment. In the densemorphology, anisotropy increases at longer length scales,while in the network morphology it is highest at shorterlength scales corresponding to narrow trails.Due to mechanical constraint from the substratum,movement of the initial bacterium is either negligible orquasi-1D due to resistance to orientational and trans-lational movement. As cell numbers increase, this 1Dconfinement facilitates collisions between cells that pro-duce enough force to change cell orientations, causingbranching. The dense configuration occurs when motil-ity cannot outpace volume expansion due to populationgrowth (Supplemental Movies S5 and S7 [5]). The diffu-sive mode corresponds to the situation where collisionsbetween small numbers of cells are sufficient to overcomethe orientational resistance applied by the substratum.In this case, branching occurs readily in arbitrary loca-tions within an existing trail and the system becomesisotropic as branches interlink (See Movie S6 [5]). Net-work morphology occurs in the intermediate case, wherecollisions between larger numbers of cells are required forbranching. The process can be observed in SupplementalMovie S4 [5].If EPS deposition is ignored and the attachment prob-ability P s is uniform, qualitatively different trail-networkpatterns are possible [Fig. S3(c,d)] that evolve from arandomly diffusing state due to the amplification of fluc- tuations in surface topography (Supplemental Movie S8[5]). These patterns are distinct from those formed byour simulations of P. aeruginosa (Fig. S4, [5]).To conclude, we have simulated biofilm growth to ex-plain how network morphogenesis in expanding bacterialcolonies arises from an interplay between mechanicalconfinement and enhanced motility due to EPS deposi-tion. These patterns cannot be explained by the EPStrail following mechanism alone, and emerge despite thetendency for exponentially growing colonies to expandisotropically. Moreover, their formation does not requirechemotaxis, contact-based signaling, surface adhesionbetween cells, or EPS-mediated nematic alignment.CZ would like to acknowledge Christian Wolff andMatthew Arnold for useful discussions. This work wasfunded by FEI Company and the Australian ResearchCouncil (DP140102721).
SUPPORTING INFORMATION FOR‘NETWORK PATTERNS IN EXPONENTIALLYGROWING 2D BIOFILMS’ by Zachreson et al.
Supplemental movies • Movie S1: The advancing edge of a
Pseudomonasaeruginosa interstitial biofilm (2000s, 0.5 frame/s). • Movie S2: The interior of a
Pseudomonas aerugi-nosa interstitial biofilm (2000s, 0.5 frame/s). • Movie S3: Simulated interstitial biofilm showingbacteria in red, EPS trails in green and substratumdeformation in gray-scale. γ = 0 . k U = 0 . . × s , 1000s/frame). • Movie S4: Simulated interstitial biofilm showingbacteria in red, EPS trails in green and substratumdeformation in gray-scale. γ = 0 . k U = 0 . . × s , 1000s/frame). • Movie S5: Simulated interstitial biofilm showingbacteria in red, EPS trails in green and substratumdeformation in gray-scale. γ = 1 . k U = 0 . . × s , 1000s/frame). • Movie S6: Simulated interstitial biofilm showingbacteria in red, EPS trails in green and substratumdeformation in gray-scale. γ = 0 . k U = 0 . . × s , 1000s/frame). • Movie S7: Simulated interstitial biofilm showingbacteria in red, EPS trails in green and substratum2deformation in gray-scale. γ = 1 . k U = 0 . . × s , 1000s/frame). • Movie S8: Simulated interstitial biofilm with EPStrail following disabled - pilus attachment is con-stant at P s = 0 . γ = 0 . k U = 0 . . × s ,1000s/frame). Anisotropy parameter calculations
For simulation data, anisotropy ∆ S φ ( r, d edge ) was cal-culated as follows: for each cell, we calculated the en-tropy of the normalized histogram of the polar coordi-nates φ defining the vectors between the cell’s centroidand the positions of k neighbors within a radius r definingthe length scale of interest as S φ ( r ) = (cid:80) i p ( φ i ) ln p ( φ i ),where i denotes one of n discrete bins in the histogram,and p ( φ i ) is the probability of randomly selecting a φ value within the i th bin.By comparing S φ ( r ) to the value associated with theuniform distribution S uni = ln n , we calculated ananisotropy parameter: ∆ S φ ( r ) = [1 − S φ ( r ) /S uni ] − ς ( n, k ) where ς ( n, k ) ≈ nk [0 . n − . + 0 . k (cid:29) n .For Figs. 5(d,e,f) and S4(a,b), the distance from thecolony edge d edge was calculated as the shortest distancefrom the cell of interest and the convex hull of the con-figuration.For analysis of experimental data, the procedure wassimilar. However, the low-resolution images of P. aerug-inosa used for comparison did not allow precise localiza-tion of cell centroids. To circumvent this limitation, theimages [such as those shown in Fig. 2(d-f)] were convertedto binary using a standard thresholding procedure fol-lowed by manual cleaning of pixels not representing bac-teria. These binary maps were then re-sampled to reflectthe number of cells (cid:104) N (cid:105) expected to occupy the coveredarea. This estimation was achieved by calculating the av-erage individual cell area a o from high-resolution imagessuch as those represented in Fig. 1(a,b), and calculating (cid:104) N (cid:105) = A tot /a o , where A tot is the total area covered bycells in the low-resolution binary image. Then, (cid:104) N (cid:105) filledpixels were randomly selected from the binary image, toapproximate a true cell configuration. To avoid artifactsdue to the square pixel grid, the position of each ran-domly selected pixel was randomized by adding a randomperturbation ∆ x ∈ [ − w/ , w/
2] where w is the averagecell width. The anisotropy of 20 such configurations wereaveraged for each of 4 replicate experiments at each gel-lan gum concentration investigated. These 4 replicateswere then averaged to produce the anisotropy plots inFig. 5(a-c). For these plots, d edge was calculated as thedistance from the cell of interest and the portion of the convex hull representing the true colony edge (ignoringthe borders of the imaging frame). To avoid edge arti-facts in the anisotropy plots, only cells further than thelength scale r from the image frame were included in theanalysis. Results of model variants lacking substratumdeformation
To substantiate the claim that EPS trail following isnot sufficient to produce the network patterns observedexperimentally, we performed tests with several modelvariants that do not include the phenomenon of substra-tum deformation. The first of these is a highly idealizedmodel that does not include cell-cell repulsion. That is,the simulated bacteria can occupy the same space. Re-sults obtained with two subvariants of this model areprovided in Fig. S1. In the first subvariant, we set themaximum pilus retraction period ( t ret = 0 . s ), 2 or-ders of magnitude faster than the experimentally derivedperiod of t ret ≈ s , used in the main text. The fast re-traction period allows the bacteria to continuously sam-ple the space in front of them, and follow EPS trails veryeffectively. Because cell-cell repulsion is disabled, the dy-namics of the bacteria are determined primarily by thedeposition rate ( k eps ) and degradation rate ( β eps ) of theEPS trails that modulate the pilus attachment probabil-ity. This is shown in Fig. S1(a), which illustrates steady-state configurations of fixed cell populations ( N = 250)in a space with periodic boundary conditions, for varied k eps and β eps values. With continuous trail sensing andno repulsive interactions between cells, trail formationoccurs over a large area of parameter space, as demon-strated by the high anisotropy values illustrated in thecorresponding interpolated contour shown in Fig. S1(b)(here, (cid:104) ∆ S φ (cid:105) r represents the spatial anisotropy averagedover length scales r ∈ [10 , , , , w for the entiresimulation duration of t f = 5 × s , measured for eachcombination of k eps and β eps ). However, if the maxi-mum pilus retraction period is increased to the experi-mentally constrained value of t ret = 10 s , trail follow-ing does not occur for any of the values of k eps and β eps tested. This is demonstrated clearly by the resulting con-figurations shown in Fig. S1(c), and the correspondinglow anisotropy values shown in Fig. S1(d).So far we have demonstrated that network formation isextremely sensitive to the trail sampling rate. While weare confident in our estimation of the true sampling ratebased on the available literature, it is possible that thisvalue is sensitive to experimental conditions. Therefore,we continue our investigating of the continuous samplingcondition in the second model variant by holding the re-traction period at t ret = 0 . s (continuous sampling), andenabling repulsive contacts between cells. Because phys-ical contact between rod-shaped cells can lead to align-3ment of orientations and collective motion, its influenceover the system’s dynamics cannot be ignored. Becausethe frequency of cell-cell interactions can be expected todepend on density, we proceed by investigating two dif-ferent cell populations.In Fig. S2(a), configurations of N = 125 cells are illus-trated and, while network formation is not as pronouncedas it is if repulsion between cells is ignored [Fig. S1(a)], itis still observable in configurations and the correspond-ing anisotropy values [Fig. S2(b)]. However, if the cellpopulation is low and nutrients are plentiful, the pop-ulation can be expected to increase. If the populationis doubled to N = 250, network formation is no longerobserved in configurations [Fig. S2(c)], and anisotropyis low for all values of k eps and β eps tested [Fig. S2(d)].It is interesting to note that while trail network forma-tion is not favored when repulsive contacts are included,the configurations in Fig. S2(c) for high values of β eps and k eps appear to demonstrate nonequilibrium clusterformation and higher degrees of collective motion thanwould be expected without EPS trails. This apparentcompetition between collective motion effects and trailnetwork formation deserves further study.To conclude this section, we have shown that trail net-work formation is not a robust consequence of EPS trialfollowing due in part to the stochastic nature of spatialsampling with T4P. Even if that limitation is relaxed,physical interactions make the phenomenon very sensi-tive to cell density due to collective motion effects. Fordilute populations of cells that can sample their environ-ment continuously, EPS deposition may be enough fortrail following phenomena. However, such behavior isnot consistent with the biological system investigated inthis work, and our model of emergent pattern formationmust include substratum deformation to explain the ex-perimentally observed phenomena. Results of model variants lacking EPS
To detail the role of EPS in colony morphogenesis, weran simulations using a model variant that lacks EPS de-position. Without EPS, pilus attachment probability P s is uniform throughout the simulation space. In this case(for fixed motility parameters as described in the modeldescription in our earlier work [14]) colony morphologydepends only on γ , k U , and P s . We tested the model’sbehavior for γ = 1 . , . , . , . k U = 0 . , . P s = 0 . , .
25. Figure S3 shows results for P s = 0 . γ = 0 . k U = 0 .
001 [Fig. S3(a) γ = 0 . P s = 0 . k U = 0 . γ = 0 .
5, without [Fig. S4(a)],and with [Fig. S4(b)] EPS trail following. The residualof these two plots [Fig. S4(c)] illustrates how the pat-terns formed due only to substratum deformation retainanisotropy over much longer length scales due to the long,tapered nature of the trails.
Image processing and growth rate estimation
The phase contrast microscopy time-series used to cal-culate growth rate was comprised of 1000 frames cap-tured at 1 frame /2 s (2000 s). Phase contrast microscopywas performed using an Olympus IX71 inverted research4microscope with a 100x 1.4 NA UPlanFLN objective,FViewII monochromatic camera and AnalySIS Researchacquisition software (Olympus Australia, Notting Hill,VIC, Australia) fitted with an environmental chamber(Solent Scientific, Segensworth, UK).The growth rates were determined by morphologicalanalysis of the individual cells using a modified versionof the BacFormatics image processing platform [19] whichwas developed in MATLAB. The calculated growth rateswere slightly faster than the actual ones due to single-pixel trimming during image segmentation. To correctfor this, we set the individual growth rate to k g = 3 . × − µ m s − in our simulations (Fig. 1c, dashed line).The source code for the beta version of the BacFormat-ics toolbox used for image segmentation and tracking inthis work is available at:https://github.com/ithreeMIF/BacFormatics CZ/.Briefly, the phase contrast images are segmented byapplying a range filter, followed by local normalization[20], low-pass frequency filtering (using a vectorized ver-sion of the butteroworth low-pass filter found here [21]),and thresholding. The resulting binary image of the cellbodies is multiplied by another binary image that is cre-ated using morphological tophat filtering to isolate animage of the spaces between cells. The resulting binaryimage effectively segments individual cells. As a finalstep, segments are excluded from this image if they donot contain local maxima from an inverted version of theoriginal phase-contrast image. This process is summa-rized schematically in Fig. S5.Cell centroid positions, orientations, lengths l , andwidths w were estimated using the regionprops() functionin the MATLAB image processing toolbox. Cell track-ing was carried out by creating a distance matrix withthe cell centroid positions in two consecutive frames andlinking cell IDs based on the minimum distance traveledbetween frames.Before estimating growth rate, we examined the lengthof each tracked cell as a function of time and identifieddiscontinuities. These were identified based on a cutofflength variation between consecutive frames ∆ l min = 0 . µ m (10 pixels) that is biologically unrealistic and is sig-nificant with respect to random variation in segment di-mensions. Such discontinuities can occur naturally whencells divide, or when there are errors in segmentation thatmerge two or more cells, or split individual cells.Between such discontinuities, the aspect ratio κ = l/ (cid:104) w (cid:105) t was calculated for each frame where (cid:104) w (cid:105) t is thetime average of the cell width, which should not changeas the cell elongates. We estimated the rate of change ofthe aspect ratio for each continuous timeseries of lengthsby applying a linear regression to each data sequence ex-ceeding a threshold duration (∆ l min = 400 s , 200 frames).The slope of the linear regression accounts for an indi-vidual growth rate value, which we averaged to find theexpectation value for growth rate in each experiment. The k g values in the main text (Fig. 1c) therefore repre-sent linear approximations of the rate at which the aspectratio of each cell is increasing. This is appropriate sinceour model uses the cell width as its fundamental unit oflength w = 1 µ m. [1] E. Ben-Jacob, O. Schochet, A. Tenenbaum, I. Cohen,A. Czirok, T. Vicsek, et al. , Nature , 46 (1994).[2] F. Farrell, O. Hallatschek, D. Marenduzzo, and B. Wa-claw, Phys. Rev. Lett. , 168101 (2013).[3] J. Howard, S. W. Grill, and J. S. Bois, Nat. Rev. Mol.Cell Bio. , 392 (2011).[4] E. S. Gloag, L. Turnbull, A. Huang, P. Vallotton,H. Wang, L. M. Nolan, L. Mililli, C. Hunt, J. Lu, S. R.Osvath, et al. , Proc. Natl. Acad. Sci. U.S.A. , 11541(2013).[5] See Supplemental Material for movies of colony morpho-genesis, details of image segmentation and anisotropy cal-culations, and a discussion of model variants lacking EPSdeposition.[6] A. B. Semmler, C. B. Whitchurch, and J. S. Mattick,Microbiology , 2863 (1999).[7] E. J. Schantz and M. A. Lauffer, Biochemistry , 658(1962).[8] V. Pelicic, Mol. Microbiol. , 827 (2008).[9] J. M. Skerker and H. C. Berg, Proc. Natl. Acad. Sci.U.S.A. , 6901 (2001).[10] B. Maier and G. C. Wong, Trends Microbiol. , 775(2015).[11] A. Gelimson, K. Zhao, C. K. Lee, W. T. Kranz, G. C. L.Wong, and R. Golestanian, Phys. Rev. Lett. , 178102(2016).[12] W. T. Kranz, A. Gelimson, K. Zhao, G. C. L. Wong, andR. Golestanian, Phys. Rev. Lett. , 038101 (2016).[13] L. Turnbull and C. B. Whitchurch, “Motility assay:Twitching motility,” in Pseudomonas Methods and Pro-tocols , edited by A. Filloux and J.-L. Ramos (SpringerNew York, New York, NY, 2014) pp. 73–86.[14] C. Zachreson, C. Wolff, C. Whitchurch, and M. Toth,arXiv preprint arXiv:1611.04223 (2016).[15] J. S. Mattick, Annu. Rev. Microbiol. , 289 (2002).[16] P. Ghosh, J. Mondal, E. Ben-Jacob, and H. Levine, Proc.Natl. Acad. Sci. U.S.A. , E2166 (2015).[17] N. M. Oliveira, K. R. Foster, and W. M. Durham, Pro-ceedings of the National Academy of Sciences , 201600760(2016).[18] D. Boyer, W. Mather, O. Mondrag´on-Palomino,S. Orozco-Fuentes, T. Danino, J. Hasty, and L. S. Tsim-ring, Phys. Biol. , 026008 (2011).[19] L. Turnbull, M. Toyofuku, A. L. Hynen, M. Kurosawa,G. Pessi, N. K. Petty, S. R. Osvath, G. C´arcamo-Oyarce,E. S. Gloag, R. Shimoni, et al. , Nat. Commun. (2016).[20] Xiong, Guanglei. (2005). Local Normal-ization. Retrieved on Aug 1, 2016 fromhttp://au.mathworks.com/matlabcentral/fileexchange/8303-local-normalization.[21] Athiq, Mohamed. (2013). Frequency domain filteringfor grayscale images. Retrieved on Aug 1, 2016 fromhttps://au.mathworks.com/matlabcentral/fileexchange/40579-frequency-domain-filtering-for-grayscale-images. β = 0.01 k p = 0.01 k p = 0.001 k p = 0.1 k p = 0.05 k p = 0.005 β = 0.001 β = 0.005 β = 0.0005 β = 0.0001 k p = 0.01 k p = 0.001 k p = 0.1 k p = 0.05 k p = 0.005 β = 0.01 β = 0.001 β = 0.005 β = 0.0005 β = 0.0001 a) c) -6-5-4-3 ln k -9 -8 -7 -6 -5 ln β -6-5-4-3 ln k -9 -8 -7 -6 -5 ln β Δ S φ r Δ S φ r b) d) β eps k eps β eps k eps FIG. S1. Results of model variants lacking substratum deformation and repulsive contacts between cells. The configurationsin (a) and the corresponding anisotropy values shown in (b) were obtained by setting the pilus retraction period t ret = 0 . s ,allowing cells to continuously sample EPS trails (illustrated in green). If trail sampling is stochastic ( t ret = 10 s ), the resultingconfigurations (c) and anisotropy (d) show a complete lack of trail network formation. k p = 0.01 k p = 0.001 k p = 0.1 k p = 0.05 k p = 0.005 β = 0.01 β = 0.001 β = 0.005 β = 0.0005 β = 0.0001 k p = 0.01 k p = 0.001 k p = 0.1 k p = 0.05 k p = 0.005 β = 0.01 β = 0.001 β = 0.005 β = 0.0005 β = 0.0001 a) c) -6-5-4-3 ln k -9 -8 -7 -6 -5 ln β Δ S φ r -6-5-4-3 ln k -9 -8 -7 -6 -5 ln β Δ S φ r b) d) k eps k eps β eps β eps FIG. S2. Results of model variants lacking substratum deformation, but including repulsive contacts between cells. Configura-tions of cells with EPS trails in green and the corresponding anisotropy values are shown for N = 125 (a, b), and N = 250 (c,d). In both cases represented, EPS trail sampling was continuous ( t ret = 0 . s ). FIG. S3. Colony morphologies simulated by the bacterial biofilm model without EPS deposition for constant surface attachmentprobability P s = 0 .
5. Surface topography after 11 replication cycles shows dense, dendritic, network, and diffusive morphologies,depending on the substratum parameters γ and k U (a, b). Dotted boxes indicate simulations where bacteria crossed the periodicboundaries before t f . Contour plots (c-f) show areas of substratum deformation where C s > C s > .
025 (d, f) forthe indicated replication cycles and substratum parameters. r d edge a) Δ S φ r d edge c) Δ S φ b) Δ S φ r d edge FIG. S4. Intensity maps of ∆ S φ as a function of length scale ( r ), and distance from the colony edge ( d edge ) for simulated biofilmswith substratum deformation only (a), EPS trail following as well as substratum deformation (b), and the corresponding residual(c), c = | a − b | . For both cases, substratum stiffness γ = 0 . k U = 0 . raw segmented local maxima range filter, local normalization, LP frequency filter inversion Morphological tophat, inversion filtered exclusion multiplication, closingraw segmented local maxima range filter, local normalization, LP frequency filter inversion Morphological tophat, inversion filtered exclusion multiplication, closing