## Abstract

Self-organization in developing living organisms relies on the capability of cells to duplicate and perform a collective motion inside the surrounding environment. Chemical and mechanical interactions coordinate such a cooperative behaviour, driving the dynamical evolution of the macroscopic system. In this work, we perform an analytical and computational analysis to study pattern formation during the spreading of an initially circular bacterial colony on a Petri dish. The continuous mathematical model addresses the growth and the chemotactic migration of the living monolayer, together with the diffusion and consumption of nutrients in the agar. The governing equations contain four dimensionless parameters, accounting for the interplay among the chemotactic response, the bacteria–substrate interaction and the experimental geometry. The spreading colony is found to be always linearly unstable to perturbations of the interface, whereas branching instability arises in finite-element numerical simulations. The typical length scales of such fingers, which align in the radial direction and later undergo further branching, are controlled by the size parameters of the problem, whereas the emergence of branching is favoured if the diffusion is dominant on the chemotaxis. The model is able to predict the experimental morphologies, confirming that compact (resp. branched) patterns arise for fast (resp. slow) expanding colonies. Such results, while providing new insights into pattern selection in bacterial colonies, may finally have important applications for designing controlled patterns.

## 1. Introduction

The collective motion rather than the migration of single individuals is known to drive the macroscopic evolution of many biological processes, such as wound healing [1,2], biofilm formation [3], tumour growth [1,4–6] and morphogenetic processes [1,7–9]. In fact, living cells and bacteria tend to form closely packed clusters (e.g. the columnar structures emerging in expanding living colonies) in which microscopic self-interactions result in cooperative migration [7,10]. Although each individual in a population can determine its own fate, recent findings highlighted that single behaviours and abilities are adjusted, thanks to stochastic differentiation, in order to fit the needs of the population as a whole [11,12].

In order to coordinate such cooperative activities, the chemical and the mechanical interactions both among individual entities and with the extracellular environment are of paramount importance. Chemical communication relies either on the capability of the cells to secrete chemicals (in cell-to-cell chemical signalling) or on the ability to sense an external chemical field, by binding specific signal molecules through their membrane receptors [13]. The chemical interaction can not only occur between adjacent cells, but can also act over long distances, through complex intracellular mechanisms involving signal transduction networks and gene network dynamics [14]. Nevertheless, living matter is not only responding to soluble biochemical signals, but also to physical factors, e.g. through surface cell receptors such as integrins and focal adhesion proteins, establishing a mechanical feedback of fundamental importance in both physiological and pathological conditions [15,16]. Thus, endogenous and exogenous physical forces act as key regulators of important intracellular signals that drive the dynamical evolution of living organisms through mechanotransduction [17].

In particular, we focus here on pattern formation in expanding bacterial clusters, studying the influence of the chemotactic motility of the colony, the interaction between the bacteria and the substrate as well as the size effects in experiments. A bacterial colony is an excellent biological system model to study general principles of collective motion and pattern formation in complex living organisms [10]. Furthermore, the study of cooperative migration and the onset of branched morphologies in bacterial colonies is an extremely multidisciplinary field of research, combining biological information with the mathematical theories of nonlinear dynamics and the physics and mechanics of non-equilibrium processes.

Thus, bacterial capability to develop elaborate branched patterns, even starting from an initially homogeneous microbial monolayer, has been intensively studied from the biological point of view, plating different bacterial colonies on a Petri dish covered with agar, under different environmental conditions [18–23]. Some of the observed patterns are reported in figure 1*a–c*, ranging from a disc-like colony (figure 1*a*) to a densely branched morphology (figure 1*c*). These biological observations highlight that nutrient diffusion and the interaction with the substrate, together with the cellular capability to proliferate and move either in a random-walk-like fashion or in response to external chemical signals (i.e. chemotaxis), are the key ingredients in the progression of front instabilities.

Indeed, both random and biased flagellation-based, run-and-tumble bacterial motions, in conjunction with collective lubrication by secretion of surfactants, enable rapid colony expansion (i.e. up to centimetres per hour [25]), as extensively studied for many bacterial species. However, we remark for sake of completeness that some bacteria also possess other motility mechanisms for colonization. For example, recent studies on *Paenibacillus dendritiformis* [25] and on myxobacteria [26] have shown that very long bacterial cells do not swarm with the standard dynamic patterns of whirls and jets, but they rather form long tracks in which each individual bacterium periodically reverses its direction, moving back and forth along moderately curved lines. The short time between this direction switching was found to be independent of the number of neighbours and the environmental factors, suggesting the existence of an extraordinary robust internal clock for reversal events in bacteria. Although the evolutionary advantages of such a reversible motion remain unclear and deserve further studies [26], we focus in the following only on the standard run-and-tumble colony expansion.

Many *in silico* mathematical models have been proposed to reproduce the spontaneous onset of the complex patterns observed during microbial growth. The theoretical approaches can be divided into two main categories: hybrid and continuous models. In hybrid models [18,27–30], the microorganisms are represented as discrete, moving entities, whereas the time evolution of the chemicals is described by reaction–diffusion equations. In continuous models [19,22,29,31–34], the bacteria as well as the nutrients and all the other possible factors involved in the process are represented via their density per unit surface. The most popular models of bacterial growth are systems of reaction–diffusion equations [22,23,31–35], allowing the description of colony patterns, including not only spreading disc-like patterns, which was found to be consistent with the solution of a two-dimensional Fisher equation for the bacterial cell density [23,35], but also branched and fractal ones [32,33]. Because bacteria expansion is modelled only by diffusion (i.e. without allowance for chemotaxis), instabilities arise either because of a nonlinear (e.g. density-dependent) diffusion coefficient [22,33] or thanks to the introduction of a bacterial transition from the active-motile and proliferative state to a passive one [31]. Furthermore, the observed branched patterns have been recovered using the assumption of a nutrition-limited process [19,22,24,33,36] or by combining signalling from a chemorepellent and a chemoattractant coupled with a proper dynamic for the colony, e.g. derived from the classical Keller–Segel model [37].

All of these studies have focused on diffusing chemicals as the major driving factors of the process, neglecting any mechanical balance laws. Some efforts to include the mechanical considerations in the modelling of bacterial expansion during the colony expansion have been done in hybrid models [28], having the major drawback of the small number of individuals that can be simulated numerically.

Moreover, the modelling of the viscous interaction between the colony and the substrate has been recently proposed in a continuous model of biofilm formation [3], proving the instability of a bacterial biofilm with planar front whose expansion relies on a nutrient-driven volumetric mass source without any chemotactic motion.

In this work, we investigate the growth and chemotactic mobility of a biological colony, coupled with the diffusion and consumption of nutrients provided by the agar, using a continuum mechanical model at the macroscopic scale. Because bacteria maintain contact with their neighbours and no gaps appear during culture expansion, the continuum assumption is effectively formulated. Moreover, bacterial growth within a Petri dish is described as a free-boundary two-dimensional problem. In the following, we first present the mathematical model based on thermomechanical considerations for describing the expansion of a circular bacterial colony. We later perform the linear stability analysis for the resulting diffusive circular growth, and we study the dynamics of pattern formation in the nonlinear regime using numerical simulations.

## 2. Mathematical model

A bacterial colony can be modelled as a two-dimensional continuum body occupying a region denoted by *Ω ^{−}*, with a moving boundary

*∂Ω*(figure 1

^{−}*d*). The bacteria are immersed into a spatial outer domain,

*Ω*

^{+}, with boundary

*∂Ω*

^{+}that stands for the border of the Petri dish with radius

*R*

_{out}. The domain

*Ω*

^{+}represents the lubricant fluid on the top of the agar, which seems to be collectively produced by the bacteria themselves, although it could also be drawn from the agar during colony expansion [18,19,24].

The mathematical model takes into account the diffusion of the nutrients in the fluid on the top of the agar, the chemotactic mobility of the bacteria and the mechanical interactions with the substrate. Here, we consider a single nutrient species with volume concentration *n*(**x**,*t*), e.g. peptone, diffusing from the outer boundary *∂Ω*^{+}, with diffusion coefficient *D _{n}* and consumed, with an uptake rate

*γ*

_{n}, only in the region occupied by the living material, such that 2.1

Considering that the living material can be macroscopically described by a Newtonian fluid moving at low Reynolds numbers [38–40], the classical Darcy's law gives the velocity of the living colony, **v** = *−K*_{p}∇*p*, where *K*_{p} is the permeability coefficient describing the friction properties with the substrate and *p* is the pressure.

Then, assuming that the colony grows, thanks to a non-convective mass flux **m**, without any significant volumetric mass source, the standard mass balance for the bacterial spatial density *ρ* reads . Being bacteria mostly composed by water, the incompressibility constraint leads to

Biological evidence suggests that bacteria are able to implement directional movements, by decreasing the tumbling frequency of their flagella, when they move up the gradient of the chemoattractant [41]. Thus, by neglecting the random motion of bacteria with respect to the directional one, we can assume that the non-convective mass flux vector is directed along the gradient of the chemical concentration, such that **m** = *χρ∇n*, where *χ* is the chemotactic coefficient [37]. Substituting Darcy's law in the mass balance equation for a homogeneous microbial colony gives the following relation between the pressure *p* and the nutrient concentration
2.2

In order to solve the system of partial differential equations (2.1)–(2.2), boundary conditions are to be provided. In particular, on *∂Ω ^{−}*, we apply the Young–Laplace equation for the pressure, the compatibility condition for the moving interface and the classical continuity conditions for the nutrient concentration and its normal derivative
2.3and
2.4where

*C*is the local curvature of the free boundary,

*σ*

_{b}is the surface tension of the interface,

*p*

_{0}is the constant outer pressure and

**n**is the outward normal vector at the boundary. It is worth noting that the surface tension of the bacterial colony results from the collective interaction with the biopolymers forming the surrounding liquid environment [42].

In the following, we apply this model to a initially circular bacterial colony of radius *R*_{0}^{*} = *R*^{*}(*t* = 0), deriving the non-dimensional system of governing equations for the dimensionless chemical concentration, and the dimensionless pressure,
2.5and
2.6using the following characteristic time *t*_{c}, length *l*_{c}, velocity *v*_{c}, pressure *p*_{c} and chemical concentration *n*_{c}: , , , , , where is the fixed concentration at the outer border of the Petri dish. The dimensionless boundary conditions for the pressure, the chemical concentration of the nutrients and the velocity, on the colony front *∂Ω ^{−}*, read
2.7
2.8
2.9On the other hand, the boundary condition on the nutrient concentration on ∂

*Ω*

^{+}is given by . For the sake of simplicity, we omit in the following the barred notation to denote dimensionless quantities.

The non-dimensionalization procedure leads to the definition of four dimensionless parameters: two of them, and are related to the chemotactic response and the bacteria–substrate interaction (*motility* parameters), whereas the other two, (i.e. dimensionless initial radius of the circular colony) and *R*_{out} (i.e. the dimensionless outer radius of the Petri dish) define the geometrical properties of the system with respect to the diffusive length, *l*_{c} (*size* parameters). In particular, the dimensionless parameter *β* represents the ratio between the energy required for the chemotactic expansion of the colony and the energy provided by the diffusive nutrients, whereas considering that *K*_{p} can be related to the friction between the colony and the substrate, *ζ*, through *K*_{p} = *l*_{c}/*ζ*, the dimensionless parameter *σ* becomes the ratio between the surface tension of the bacterial colony and the product of the colony–substrate friction and the diffusion coefficient, i.e. *σ* = *σ*_{b}/(*D*_{n}*ζ*).

## 3. Results

### 3.1. Linear stability analysis

The linear stability analysis for standard problems with a moving interface is usually performed by perturbing a steady state or a steady wave profile in a comoving frame, such as travelling wave solutions for infinite rectilinear fronts. Because the governing equations for the finite bacterial domain do not admit such kinds of solution for a circular front, let us consider the quasi-stationary expansion of a circular bacterial colony of radius *R*^{*}(*t*), assuming that the diffusive process is much faster than the border expansion, so that the time-derivative term in (2.5) can be neglected. The quasi-stationary velocity of the circular front has only a radial component, i.e. , with , where *I _{j}*(

*r*) is the modified Bessel function of the first kind of order

*j*, evaluated in

*r*, and is the nutrient concentration at the free boundary, in the quasi-stationary state. Further details on the derivation of the quasi-stationary solution are presented in the electronic supplementary material.

Considering a perturbation of the moving interface given by *R*(*θ*, *t*) = *R**(*t*) + *ɛ*e^{λt+}^{i}* ^{kθ}*, it is possible to find the dispersion equation relating the time growth rate

*λ*to the integer wavenumber

*k*in an implicit way, as a function of the four dimensionless parameters

*β*,

*σ*,

*R** and

*R*

_{out}3.1for

*λ ≠ −*1. More information on the procedure used to obtain the dispersion relations and the determination of the coefficient

*A*

_{λ}can be found in the electronic supplementary material.

Interestingly, equation (3.1) shows that the amplification rate of the perturbation over time depends on the radius of the colony, *R**. The dispersion diagrams presented in figure 2, obtained solving iteratively the dispersion equation (3.1) through the secant method, demonstrate that the colony boundary is always linearly unstable for large wavelengths. The resulting dispersion curves show a stabilization of the front at short-wavelengths, which is also typically found in some classical fluid instabilities [43]. Indeed, the governing equations closely resemble those driving the Saffman–Taylor instability in a Hele–Shaw cell, where *σ* has the same physical meaning of the capillary number, introducing a correction in the dispersion relation proportional to (*k* − *k*^{3}) [44]. Such a regularizing effect is even enforced in our model by the presence of the Laplacian operator at the right-hand side of equation (2.2), corresponding to a mass source term which is not considered in viscous instabilities. While for a planar front, one would expect to find *λ* = 0 for *k* = 0 for translational symmetry, we note that for the circular front we consider , because *k* = 0 would correspond to the quasi-static expansion of the circular front. Furthermore, even if *β* somewhat physically corresponds to the influx volume flow rate in an initially circular front, the velocity of the quasi-stationary front depends in this problem on the initial radius *R**, while it is often taken as a constant in viscous fingering problems [45], thus introducing a marked size-dependence in the resulting dispersion equation. In both cases, the linear stability analysis can identify the finite critical (or characteristic) wavelength with the most unstable growth rate, which increases for increasing *β*. Nonetheless, even if such linear stability curves are indicative of unstable front dynamics, suggesting the nonlinear formation of dendritic patterns [46], we also expect to find differences in the fully nonlinear dynamics owing to the more complex functional dependence on the growth and geometrical parameters in this problem. Furthermore, the maximum time growth rate increases as far as *β* increases, as shown in figure 2*b*,*d*. Comparing the dispersion curves for different dimensionless radii of the colony, while maintaining the ratio between *R** and *R*_{out} as a constant (which corresponds to a variation of the diffusive length), it is possible to note that increasing such a size parameter increases the range of unstable wavenumbers and the characteristic wavenumber, but it decreases the value of the fastest growth rate (figure 2*b* versus *d*).

### 3.2. Numerical simulations

Because the linear analysis indicates the occurrence of a front instability with a finite characteristic wavenumber, we investigate the nonlinear pattern formation of the bacterial colony using a finite-element code implemented in FreeFem++ (http://www.freefem.org). The equations (2.5)–(2.6) are solved on an adaptively refined triangular grid, fitting the moving interface. Given the concentration of nutrients at time *t _{i}*, we first compute the pressure

*p*, through equation (2.6) and then the velocity field, using Darcy's law. We are yet able to explicitly move the boundary and solve (2.5) for the concentration at time

_{i}*t*

_{i}_{+ 1}, using an implicit-Euler scheme.

The numerical results are reported in figure 3*a* showing the emerging patterns versus the two dimensionless parameters *β* and *R*_{out}, at fixed and *σ*.

In particular, it is found that the branching of the colony is strongly favoured by smaller values of the motility parameter *β*. Because *β* drives the velocity of the quasi-stationary front, we find that the instability of the contour is enhanced by a low expanding velocity, which is a common feature of growing biological systems in a diffusion-dominated regime [5]. In fact, small values of *β* correspond to situations in which the diffusion of chemicals is dominant on the chemotactic mobility of the colony.

Furthermore, figure 3*a* points out the existence of a size effect in pattern formation, showing that the size parameter *R*_{out} determines the characteristic wavelength of the branching process. In fact, keeping the aspect ratio constant, a larger Petri dish *R*_{out}, with respect to the diffusive length *l*_{c}, determines branched patterns characterized by a smaller wavenumber.

The morphological diagram in figure 3*a* also reveals that the number of emerging fingers perfectly corresponds to the critical wavenumbers predicted by the linear stability analysis. Indeed, for the particular case *β* = 1, and *R*_{out} = 155, the blue dispersion curves in figure 2*a,b* predict the fastest linearly unstable mode with *k*_{c} = 10, whereas in the numerical simulations (figure 3*a* top-right), 10 branches can be observed. Moreover, for *β* = 1, and *R*_{out} = 350, the analytical characteristic wavenumber (see the blue dispersion curves in figure 2*c,d*) is *k*_{c} = 19, which is confirmed by the numerical simulation in figure 3*a* bottom-right, where 19 totally developed branches appear. Considering *β* = 4.25, the maximum in the red dispersion curve in figure 2*d* shifts to *k*_{c} = 22, which is in agreement with the 22 fingers depicted in figure 3*a* bottom-centre, where the same parameters are considered. The agreement between numerical and analytical results is good whenever the velocity of the front is slow enough to satisfy the quasi-stationary hypothesis, suggesting that in this regime the nonlinear effects fix the linearly unstable pattern. Conversely, the nonlinearities have a much stronger stabilizing effect on the colony profile when considering initially fast front dynamics, corresponding to high values of the parameter *β*. This fully nonlinear behaviour is strikingly similar to the diffusion-dominated instability in other biological systems, where rounder patterns are observed for increasingly faster initial fronts [5].

The morphological diagram in figure 3*a* also highlights that the occurrence of branching can be detected when the ratio between the area of the colony and its perimeter over time (square markers) deviates from half of the average radius of the colony (dashed lines).

Another important parameter in the definition of branched patterns is the roughness of the profile, measured as the root of the mean square deviation of points on the front from the average radius of the colony [29]. Either plotting this parameter as a function of time (figure 3*a*) or as a function of the averaged radius of the colony, , normalized with respect to the initial radius, (figure 3*b*), it is possible to observe that the interface roughness continuously increases for branched patterns, whereas it later saturates to an almost constant value only in compact colonies, in agreement with the experimental measurements performed in reference [47]. In summary, it is found that the motility parameter *β* defines the occurrence of the branching regime, whereas the size parameter *R*_{out} determines the characteristic wavelength of the branching process.

Furthermore, numerical simulations agree with the theory in predicting that an increase of the surface tension *σ* has a stabilizing effect also in the nonlinear regime (figure 4*a* and the electronic supplementary material, movies S3–S5). Another important parameter in the model is the initial size of the colony with respect to the diffusive length. In order to assess the influence of this parameter on the onset of branches, we have run a set of simulations keeping *R*_{out}, *σ* and *β* fixed and letting vary. The resulting morphological diagram is reported in figure 4*b*, showing that more fingers develop as the initial radius of the colony increases.

Let us now briefly discuss the characteristics of the branching process in the expanding colony. Once branching has been triggered, all developing fingers elongate in the radial direction. Further branching can occur, even though some of the second-generation fingers remain very short because of the geometrical constraints imposed by their neighbours. In figure 5, the geometrical characteristics of the fingers are plotted versus time. We define the finger base, *B*, as the distance between two subsequent points of local maximum for the curvature of the moving front. The amplitude of the finger, *A*, is defined as the maximum distance between the points belonging to the finger contour and the corresponding finger base segment. The numerical results in figure 5 correspond to the simulation on the upper left of figure 3. It is observed that at an early stage (up to *t* = 1650) both the amplitude and the base strongly increase for each of the five evolving fingers. As soon as second-generation fingers appear (from *t* = 1650 on), the amplitude of the new ones strongly increases, whereas the base remains almost constant in time. We report in the table within figure 5 the parameters of the power-law curve, *c · t ^{α}*, which best fits the numerical data for the base and the amplitude of the fingers. Interestingly, at the early stage, the best-fitting exponent for the ratio amplitude/base of the fingers is ≈0.45, which is very close to the expected square root growth over time exponent in a diffusion-driven instability [43].

Let us now return to the dimensional physical quantities in order to discuss the results of the numerical simulations with respect to the biological data. Available data in experimental literature for expanding colonies are shown in table 1. In particular, in the following, we compare the results of the simulated patterns of our model, reported in figure 1*d–f*, with their corresponding experimental counterparts in figure 1*a–c*. In all cases, the Petri dishes in the experiments have the standard radius of 44 mm.

As shown in figure 1*d*, the evolution of the disc-like colony observed in reference [24] can be reproduced by taking a diffusion coefficient *D _{n}* = 3 × 10

^{−}^{11}m

^{2}s

^{−1}, and a small uptake rate,

*γ*= 3.7 × 10

^{−}^{4}s

^{−}^{1}, giving a characteristic length

*l*

_{c}≈ 285 μm, corresponding to a dimensionless external radius of

*R*

_{out}≈ 155. Fixing

*β*= 8.5, we obtain an average velocity of the front equal to ≈3.2 mm h

^{−1}, which is in the order of magnitude of the values found in literature for round and compact colonies, such as the one in figure 1

*a*(table 1). Considering the characteristic concentration

*c*= 0.1 mM, as the one used in reference [48], the parameter

_{n}*β*= 8.5 describes a chemotactic coefficient

*χ*≈ 48

*×*10

^{−}^{5}cm

^{2}/(s mM), which is also in the biological range (table 1).

Thus, disc-like patterns are found to arise in situations in which the chemicals have an intermediate diffusion rate and are consumed by bacteria at a slow rate. Moreover, an Eden-like pattern (i.e. a clusters whose inner structure is almost completely compact but whose surface is comparatively rough [23]) is observed considering a smaller diffusion coefficient *D _{n}* = 10

^{−}^{11}m

^{2}s

^{−1}[24] and an uptake rate

*γ*= 6.5 × 10

^{−}^{4}s

^{−}^{1}(i.e.

*l*

_{c}≈ 124 μm and

*v*

_{c}≈ 290 μm h

^{−1},

*R*

_{out}≈ 350) as depicted in figure 1

*e*. Keeping

*β*= 8.5, the resulting mean front velocity is equal to 136 μm h

^{−1}, which is in agreement with the reported velocities for such patterns (table 1). In this case, under the same assumptions used above, we obtain a consistent chemotactic coefficient of about

*χ*≈ 0.85

*×*10

^{−}^{5}cm

^{2}/(s mM). Therefore, our model predicts that Eden-like structures can be obtained in the range of relatively small diffusion coefficients of the nutrients, intermediate uptake rates and small chemotactic behaviour. Although reproducing the typical macroscopic roughness of a compact colony, we have to remark that our continuous model is not able to capture the observed microscopic roughness [31], which would require a discrete approach at a cellular scale.

Furthermore, we have found that branched structures arise in simulations with high values of the dimensionless external radius *R*_{out} and small values of the motility parameter *β*, physically corresponding to high values of both nutrient diffusion and uptake rate, coupled with a low concentration of nutrients in the agar. Indeed, the development and evolution of the fingers reported in reference [13] are reproduced by our simulations in figure 1*f*, considering *γ* = 0.11 s^{−}^{1} and *D _{n}* = 10

^{−}^{10}m

^{2}s

^{−1}, so that

*l*

_{c}≈ 30 μm and

*t*

_{c}≈ 9 s. Accordingly, the observed displacement of about 2.4 mm in the first 22 h, for an initial colony with a diameter of about 6 mm [13], perfectly corresponds to the displacement of the front recorded in this simulation (i.e. 82

*·l*

_{c}at

*T*= 9000, for a colony with

*R*

_{0}* = 100). In this case, we fix

*β*= 1, corresponding to a chemotactic coefficient

*χ*≈ 10

^{−}^{5}cm

^{2}/(s mM), which is slightly below the typical values reported in reference [48].

In summary, our model is able to predict the observed experimental morphologies, reproducing that higher velocities of the front, which are linked to a stronger response of bacteria to the external chemical field (i.e. high *χ*), correspond to more rounded profiles, in accordance with biological observations [22,31].

Let us finally add a few considerations about the role of surface tension for pattern formation. Considering a friction coefficient *ζ* = 1 nNs/(μm) (compatible with the values reported in table 1), the surface tension for the three simulations results in the range *σ*_{b} ≈ 0.09−0.7 nN μm^{−1}, which is slightly higher than the values found for cell clusters (table 1). Accordingly, although we find that disc-like colonies can also be numerically obtained by increasing the dimensionless motility parameter *σ* (figure 4), this would correspond to much higher values of the surface tension than the biological range. In conclusion, we argue that the compact expansion of the bacterial front likely relies more on the capability of bacteria to actively move and consume nutrients and on the diffusive properties of the chemicals in the agar, rather than on the surface tension and the adhesive properties of the colony.

## 4. Discussion

In this work, we have presented an analytical and computational analysis of a continuum model for studying pattern formation during the spreading of an initially circular bacterial colony on a Petri dish. The proposed model differs from previous ones [3,23,33] by taking into account both the chemical effects, dictated by the diffusion of nutrients and the chemotactic response of bacteria, and the mechanical/viscous interaction between the colony and the substrate. In particular, four dimensionless parameters are found to characterize the model dynamics: two of them describe the factors that drive the *motility* of the colony, i.e. *β* represents the competition between chemotactic and diffusive effects, whereas *σ* is the ratio between the surface tension of the colony and the friction with the substrate, whereas the initial dimensionless radius of the colony, *R*_{0}*, and the dimensionless radius of the Petri dish, *R*_{out}, account for *size* effects.

The linear stability analysis has evidenced that an initially circular colony is always linearly unstable to perturbations of the interface, having the typical dispersion curves found for fluid instabilities, such as viscous fingering. The numerical simulations of bacterial spreading in the fully nonlinear regime have confirmed the development of undulated finger-like structures, which align in the radial direction and later undergo further branching, and the existence of a characteristic wavenumber, as predicted by the linear stability analysis. While the finger dynamical behaviour shows the occurrence of the typical shielding and tip-splitting found in Saffman–Taylor instability [44], we find here a more complex dependence of the characteristic length scales of such fingers on the size parameters of the problem, with more branched patterns appearing with larger Petri dishes, whereas the emergence of branching is driven by the motility parameters, being especially favoured by small values of *σ* and *β*, which also correspond to slow moving fronts. However, through comparison of the model parameters with real biological data, we argue that the selection mechanism of branched patterns likely relies on the *motility* parameter *β*, i.e. on the interplay between diffusion and chemotaxis, rather than on *σ*. In particular, in such a range, the aspect ratio of the fingers is found to scale to approximately the square root of the time. Because *β* drives the velocity of the colony interface, the model also confirms the experimental observations that compact (resp. branched) patterns arise for fast (resp. slow) expanding colonies [22,31]. Furthermore, the scaling laws for the roughness coefficient of the interface in numerical simulations are consistent with the ones reported in the biological investigations [35].

Interestingly, the modelling of the branching structures does not require neither the introduction of a nonlinear (e.g. density-dependent) diffusion coefficient, as in references [22,33,54], nor the transition from the active-motile and proliferative state of bacteria to a passive state, as in reference [31].

The experimental morphologies reported in references [13,23,24] (figure 1*a–c*) are also compared with the patterns resulting from the numerical simulations of the proposed model (figure 1*d–f*), where the relevant biological parameters in each experimental setting are used as input. A striking similarity has been found between experiments and model simulations either on the emerging disc-like, Eden-like and branched structures, or on the pattern dynamics during the expansion process.

It is worth noting that the proposed model is not able to reproduce the evolution of colonies in which fractal patterns or microscopic roughness appear, because other effects should be taken in account at scales where a continuous representation of the colony is no longer a valid assumption. Alternatively, highly fractal patterns can be reproduced by using either reaction–diffusion models with nonlinear terms [24,32,33] or discrete/hybrid models [18,29].

In conclusion, the results of this study suggest that the pattern selection and evolution in expanding bacterial colonies is driven by a complex interplay of the chemotactic response, the substrate–bacteria interaction and the size effects. Furthermore, the present study highlights the necessity to perform new experimental tests with a better quantitative characterization of both the chemical response and the mechanical forces involved during the bacterial expansion, as done for eukaryotic cells [55,56]. Indeed, even if the formation of spatial patterns in growing bacterial colonies has been intensely studied during the last years, further progress in the field has been hindered by the lack of detailed experimental data [10]. Such a combined approach has the potential to foster our understanding of pattern selection and dynamics in bacterial colonies, with important applications for designing and engineering controlled patterns.

Future works will be focused on the refinement of the proposed model, simulating more realistic initial and boundary conditions, possibly supported by novel biological experiments and coupled with a more complex representation for the living material (e.g. including viscoelasticity [57]). Further developments should also include the active motility and the chemomechanical interactions occurring at the microscopic scale inside the colony, which are of fundamental importance in the appearance and evolution of patterns in some bacterial strains [10,25,26].

Finally, future applications of this model shall concern free-boundary problems in other biological systems (e.g. wound healing), aiming at giving insights on the role played by physical forces in directing self-organization during the development of living organisms.

## Funding statement

This work was partially supported by the ‘Start-up Packages and PhD Programme’ project, co-funded by Regione Lombardia through the ‘Fondo per lo sviluppo e la coesione 2007–2013-formerly FAS' and by the ‘Progetto Giovani GNFM 2014′, funded by the National Group of Mathematical Physics (GNFM-INdAM).

## Acknowledgement

We are grateful to Davide Ambrosi for helpful discussions.

- Received November 21, 2014.
- Accepted January 13, 2015.

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