The effects of nutrient chemotaxis on bacterial aggregation patterns with non-linear degenerate cross diffusion
TThe effects of nutrient chemotaxis on bacterialaggregation patterns with non-linear degenerate crossdiffusion
J. Francisco Leyva a , Carlos M´alaga b , Ram´on G. Plaza c, ∗ a Posgrado en Ciencias Matem´aticas, Facultad de Ciencias, Universidad Nacional Aut´onomade M´exico, Circuito Exterior s/n, Ciudad Universitaria, C.P. 04510 M´exico D.F. (Mexico) b Departamento de F´ısica, Facultad de Ciencias, Universidad Nacional Aut´onoma deM´exico, Circuito Exterior s/n, Ciudad Universitaria, C.P. 04510 M´exico D.F. (Mexico) c Departamento de Matem´aticas y Mec´anica, Instituto de Investigaciones en Matem´aticasAplicadas y en Sistemas, Universidad Nacional Aut´onoma de M´exico, Apdo. Postal 20-726,C.P. 01000 M´exico D.F. (Mexico)
Abstract
This paper introduces a reaction-diffusion-chemotaxis model for bacterial aggre-gation patterns on the surface of thin agar plates. It is based on the non-lineardegenerate cross diffusion model proposed by Kawasaki et al. [1] and it includesa suitable nutrient chemotactic term compatible with such type of diffusion.High resolution numerical simulations using Graphic Processing Units (GPUs)of the new model are presented, showing that the chemotactic term enhances thevelocity of propagation of the colony envelope for dense-branching morphologies.In addition, the chemotaxis seems to stabilize the formation of branches in thesoft-agar, low-nutrient regime. An asymptotic estimation predicts the growthvelocity of the colony envelope as a function of both the nutrient concentrationand the chemotactic sensitivity. For fixed nutrient concentrations, the growthvelocity is an increasing function of the chemotactic sensitivity.
Keywords:
Nutrient chemotaxis, Bacillus subtilis, nonlinear degeneratediffusion, front propagation.
1. Introduction
The spontaneous emergence of complex patterns is ubiquitous in nature.When faced against hostile environmental conditions, for example, bacterialcolonies may exhibit very diverse morphological aggregation patterns. Suchconditions can be reproduced during in vitro experiments on the surface of thin ∗ Corresponding author. Tel.: +52 (55) 5622-3567. Fax: +52 (55) 5622-3564.
Email addresses: [email protected] (J. Francisco Leyva), [email protected] (Carlos M´alaga), [email protected] (Ram´on G. Plaza)
Preprint submitted to Elsevier November 9, 2018 a r X i v : . [ q - b i o . CB ] J a n igure 1: Depiction of the morphological diagram of the aggregation patterns exhibited bycolonies of B. subtilis as a function of the concentration nutrient C n (vertical axis) and thesolidity of the agar medium expressed as 1 /C a (horizontal axis). Both concentrations aremeasured in grams per liter [g/l]. This diagram was redrawn from the one in [1]. agar plates by using a very low level of nutrients, a hard surface (high concen-tration of agar), or both. Some species (like the bacterium Bacillus subtilis )are known to exhibit a great variety of branching patterns depending on thenutrient concentration, softness of the medium, and temperature. For example,when inoculated on a nutrient-poor solid agar, the patterns show fractal mor-phogenesis similar to diffusion-limited aggregation (DLA) (see, e.g., Matsuyamaand Matsushita [2], and Ben-Jacob et al. [3]). When inoculated on semi-solidagar (softer medium) the bacterial colonies exhibit a dense-branching morphol-ogy (DBM) with a smooth colony envelope propagating outward (cf. Ohgiwari et al. [4]). Finally, when both nutrient and agar’s softness further increase,the bacteria form simple circular patterns growing homogeneously (cf. Wakita et al. [5]). A detailed account and comparison of the different experimentalresults can be found in Kawasaki et al. [1] (see also Golding et al. [6] and thereferences therein). These experimental observations are usually summarized ina morphological diagram like the one shown in Figure 1. The horizontal axismeasures the medium softness expressed as 1 /C a , where C a is the agar con-centration, and the vertical axis is the concentration of the nutrient C n . Whenthe medium is very hard, the patterns are DLA-like (region A) at low nutrientconcentration, and round and compact with fractal boundary for high values of C n (region B). DBM-like patterns arise in regions where the nutrient is poor2nd the agar softness is intermediate (semi-solid), with an envelope smoothlyrounded and moving outward (region E). In region D, where the agar is softerand the nutrient is more abundant the dense branches fuse together forminghomogeneous colonies with smooth boundaries. Although susceptible of furtherimprovements (see, for example, Wu, Zou and Jin [7]), these diagrams provide agood macroscopic description of the colony patterns in terms of simple physicalincubation conditions.In order to explain the experiments, many mathematical models have beenproposed. One method of modeling is the continuous deterministic reaction-diffusion approach in which the bacterial density and the nutrient concentra-tion are described with continuous time evolution systems of partial differen-tial equations (cf. Murray [8]). Motivated by the experimental observations ofOhgiwara et al. [4] for DBM-like transition patterns between regions E and D inthe morphology diagram (see section 2.1 below), Kawasaki et al. [1] proposedin the mid-90’s a reaction-diffusion model with non-linear density-dependentdegenerate (cross) diffusion coefficient, which captures many of the pattern fea-tures found experimentally. Thanks to its capability of reproducing the complexdense morphology of patterns as well as its rich mathematical structure, thismodel has established itself as a well-known non-linear diffusion system in theliterature (see, e.g., [8, 9, 10] and the references therein).There are, however, other phenomena which must be taken into consider-ation. Bacterial chemotaxis (that is, the process by which bacteria migratetoward higher concentrations of certain chemical fields) has attracted signifi-cant interest due to its critical role in pattern formation (see [11, 12, 13, 14] and[15] for a review). Ben-Jacob and co-workers (cf. [16, 6, 17]) have suggestedthat the chemotactic response of the bacteria is essential to understand someof the fundamental features of the observed patterns. In particular, Golding etal. [6] discussed the experimental evidence and provided theoretical argumentsfavouring the introduction of chemotaxis towards high concentration of nutri-ents into many of the the reaction-diffusion models existing in the literature.Although the authors did not propose a new chemotactic term for the particularmodel of Kawasaki and collaborators, they certainly discussed what the relationbetween a generic diffusion term and the appropriate chemotaxis should be (seesection 2.2 below).Following these ideas, this paper incorporates a suitable nutrient chemotac-tic term into the original model by Kawasaki et al. , and explores the effectsthat such chemical signaling brings upon the aggregation patterns. For thatpurpose, we have conducted high resolution numerical simulations of the newsystem of equations using Graphic Processing Units (GPUs) to perform par-allel computations. Moreover, we were able to reproduce the results withoutchemotaxis and to compare them to the new emerging patterns when the chem-ical signaling is switched on. One of the main features of the new chemotacticterm is an increase of the front propagation velocity for the colony envelope.Such observation is justified using detailed front asymptotic calculations of thespeed and the numerical results. We find that, except for a regime of soft agarand low nutrient, the chemotaxis does not modify the pattern colony structure3ignificantly, providing evidence of the conjecture of Cohen et al. [16].The remainder of this article is organized as follows. Section 2 gives ac-count of the original non-linear reaction-diffusion (non-chemotactic) model andintroduces the chemotactic term into the equations. The resulting system isfurther normalized in order to reduce the number of parameters. Section 3 con-tains the results of our numerical simulations, as well as a comparison of theobserved patterns to those without chemotaxis. In order to explain the effectsof the chemotactic terms, section 4 contains an approximation of the speed ofthe envelope front using the ideas of geometric front propagation. Furthermore,we provide a numerical estimation of the speed which approximates well thetheoretical result, and compare it to the speed without chemotaxis. The es-timations are supported by some theoretical calculations which can be foundin Appendix A at the end of the paper. Finally, in section 5, we make someconcluding remarks.
2. Mathematical modeling
In this section we introduce the mathematical reaction-diffusion-chemotaxismodel studied in this paper. We start by providing a detailed explanation ofthe original model in [1].
The model proposed by Kawasaki et al. [1] is a reaction-diffusion continu-ous system of partial differential equations that takes into account the detailedmovement of the bacteria. The model has the following general form: n t = D n ∆ n − f ( n, b ) (1a) b t = ∇ · ( D b ∇ b ) + θf ( n, b ) (1b)where n = n ( x, y, t ) and b = b ( x, y, t ) represent the concentration of the nutrientand the density of the bacterial cells, respectively. We will denote every pointin the domain as ( x, y ) ∈ Ω ⊂ R . Here Ω is a bounded domain and t ≥ D b and D n , respectively. The coefficient D n > D b , which depends on both b and n . Its design was motivated by thework of Ohgiwari et al. [4] who investigated the patterns formed by B. subtilis on thin agar plates. They discovered that variations in the environmental con-ditions as well as in the concentrations of both nutrient and agar medium leadto drastic changes in the morphology of the colony patterns. They found twodistinct type of growing processes. One is performed by immobile cells whenthe growth occurs on hard agar plates. In this case, the outermost part of thecolonies grows by cell division and cells in the inner part change into a dormantstate. The other type of growth takes place for intermediate and low values of C a . This means that the growth is produced by the movement of bacterial cells4nd, due to the softness of the agar, they can move around dynamically. Theoutermost front of the colonies is a boundary layer of cells whose movementis dull. But cells inside the layer move around actively and sometimes theybreak through the layer and immediately become dull. The bacteria inside thecolonies are usually inactive.In view of such experimental observations, Kawasaki and co-workers conjec-tured that the bacteria are immotile when either n or b are low and becomeactive as n or b increase. Therefore, they proposed D b as a non-linear (cross)diffusion function of the form D b = σnb, (2)where σ = σ (1 + ∆). The parameter σ > f ( n, b ) is the consumption rate of the nutrient by thebacteria and θf ( n, b ) is the growth rate of cells, being θ > f is of Michaelis-Menten type, f ( n, b ) = κnb γn , (3)where κ > γ > f ( n, b ) = κnb. (4)The latter is the form of the kinetic function used by Kawasaki and collaboratorsin their simulations, and the function that we consider in this paper as well. Golding et al. [6] argue that the growth of bacterial colonies is an auto or-ganization process which involves several mechanisms of chemotactic responses.They suggested that there are at least three different types of chemotactic sig-nals. One of the signals is the result of nutrient (or food) attractive chemotaxisand it is the dominant signal in certain regimes of the nutrient level. Anothertype of signal is a long range repulsive one, which is produced by the starvingbacteria in the center of the colonies. The last type of signal is secreted by thebacteria in the front that are immersed in waste products. These bacteria askfor help to metabolize waste products by segregating a chemoattractant. Therange of this signal is determined by the diffusion coefficient and the rate ofdecomposition. It is, in general, of short range.When refering to bacterial colonies, Golding et al. [6] suggested that nutrientchemotaxis is the mechanism responsible of the increment in the growth speedand the control (or even the decrease) of the fractal dimension of the arisingpatterns. Such process should permit to increase the propagation of the front.5n view of these observations, we introduce a nutrient chemotactic term intosystem (1). For that purpose, we need to define the chemotactic flux J c which,in general, is written as J c = ξ ( b ) χ ( n ) ∇ n. (5)Here χ = χ ( n ) is the chemotactic sensitivity function and ξ = ξ ( b ) is the bac-terial response to the nutrient gradient. The chemotactic signals are detectedby dedicated membrane receptors which perform a binding between ligands (at-tractors) and specific binding proteins [18]. Bacteria sense temporal gradients,that is, they compare sequential registers of their occupied chemoreceptors [19].This complexity is built into the model equations through a signal-dependentchemotactic sensitivity function. One of the most commonly used forms is the“receptor law” proposed by Lapidus and Schiller [20], χ ( n ) = χ K d ( K d + n ) (6)where χ > K d > K d for the attractor-receptor molecular interaction has a unique value that dependson the bacterial strain, nutrient type and other in vitro conditions, and has tobe determined experimentally (see Goldman and Ordal [26] for values of K d on various substrates of B. Subtilis ). Moreover, it can be interpreted as thenutrient level where the chemotactic sensitivity χ ( n ) is maximum. Thus, weshall use K d as the normalization factor for the nutrient density n .The bacterial response function ξ = ξ ( b ) is positive for attractive chemotaxisand negative for a repulsive one. In [6] the authors suggest that, if the bacteriamove within a liquid and the colony density is low, then the bacterial responseshould be proportional to the bacteria concentration times the diffusion, thatis, | ξ ( b ) | ∝ bD b . (7)Notice that, on account of (7), if the diffusion coefficient is constant then thebacterial response function is proportional to b , recovering in this case the clas-sical Keller-Segel chemotactic flux [11, 12]. Therefore, the authors suggest that,for general non-linear diffusion functions D b , the expression of the chemotacticflux J c must be modified according to (7) in order to incorporate the densitydependence of the bacterial movement. For example, for the non-linear diffusion6onsidered by Kitsunezaki [27], namely D b = D b m where D is constant and m > ξ ( b ) = bD b = D b m +1 as the appropriateresponse function. We now introduce a suitable chemotactic term into model (1) expressing thechemical signaling of the nutrient concentration n . In view of the form of thenon-linear diffusion (2) and taking into account (7), we propose the following(attractive) chemotactic flux J c = σnb χ ( n ) ∇ n, (8)where χ ( n ) is the receptor law (6). Substracting the divergence of such flux tothe bacterial density equation in (1) we obtain the following reaction-diffusion-chemotaxis model: n t = D n ∆ n − κnb (9a) b t = ∇ · ( σnb ∇ b ) + θκnb − ∇ · (cid:18) σnb χ K d ( K d + n ) ∇ n (cid:19) , (9b)for ( x, y ) ∈ Ω, t ≥
0. The system is subject to initial conditions modeling anuniformly distributed initial constant concentration of nutrient and the initialinoculum of the bacteria. They have the form n ( x, y,
0) = ˆ n , b ( x, y,
0) = ˆ b ( x, y ) . (10)Here ˆ n is a positive constant and ˆ b is a given function. The system (9) isfurther endowed with no-flux boundary conditions, ∇ b · ν = 0 , ∇ n · ν = 0 , ( x, y ) ∈ ∂ Ω , (11)where ν ∈ R , | ν | = 1, is the outer unit normal vector at the boundary of thedomain.System of equations (9) adds the nutrient chemotactic signaling to model (1).The flux (8) can be interpreted as the cross chemotactic function correspondingto the nonlinear diffusion (2). In order to reduce the number of parameters in system (9) we introduce thefollowing rescaling of the variables:˜ x = (cid:16) θK d κD n (cid:17) / x, ˜ y = (cid:16) θK d κD n (cid:17) / y, ˜ t = ( θK d κ ) t, ˜ n = nK d , ˜ b = bθK d , ˜ σ = (cid:16) θK d D n (cid:17) σ. n t = ∆ n − nb (12a) b t = ∇ · ( σnb ∇ b ) + nb − χ ∇ · (cid:18) σnb (1 + n ) ∇ n (cid:19) , (12b)with σ = σ (1 + ∆), and with initial conditions n ( x, y,
0) = ˆ n /K d ≡ n ,b ( x, y,
0) = ˆ b ( x, y, / ( θK d ) ≡ b ( x, y, , (13)for ( x, y ) ∈ Ω. To complete the system we impose no-flux boundary conditionson the rescaled variables ∇ b · ν = 0 , ∇ n · ν = 0 , ( x, y ) ∈ ∂ Ω . (14)The only remaining free parameters are σ >
0, measuring the hardness ofthe agar medium (large σ for low agar concentration); n >
0, measuring therelative initial nutrient concentration; and χ ≥
0, which measures the intensityof the chemotaxis.In the sequel we compute numerically the solutions to system (12) underconditions (13) and (14). We are particularly interested in the interplay betweenthe different parameter values expressing agar concentration, initial nutrient andthe chemotactic signal, and how they affect the colony patterns.
3. Numerical simulations
We performed several numerical simulations of the system (12) using anexplicit second order Runge-Kutta and finite differences numerical scheme andthe NVIDIA CUDA libraries for the C++ programming language. The librariesare written to run on GPUs, and the CUDA programming language is designedunder the simple program/multiple data (SPMD) programming model (cf. [28]).More specifically, the GPU architecture is very well-suited to address problemsthat can be expressed as data-parallel computations, that is, the same program(finite difference calculations) is executed on many data elements (grid points)simultaneously (see [29]). Hence, we overcome the inconvenience of small timesteps required by explicit numerical schemes for solving stiff equations, and wereable to accomplish thousands of iterations in a couple of hours. Our numericalexperiments were performed on a NVIDIA Tesla c (cid:13) C2070 graphics card with448 CUDA cores.The calculations were performed on a two dimensional square domain ofside L = 680. The selected mesh for all our computations is a square grid of N = 2048 points per side (see discussion on section 3.2 below). Therefore, thegrid width that we used was ∆ x = ∆ y = L/N ≈ . t = 0 . t = 0 . n σ χ
0, 2.5, 5.0, 7.5, 10.0
Table 1: Free parameter values used during the numerical simulations.
Following [1] and for the sake of comparison, the initial distribution for thebacteria was taken as b ( x, y,
0) = b M e − ( x + y ) / . (15)where b M = 0 .
71 is the maximum density at the center of the domain. Recallthat the initial condition for the nutrient is a fixed constant, for which we takethe values n = 0 .
71 and 1 . σ is perturbed fromthe mean σ through a random variable ∆ taken from a triangular distributionwith support on [ − , σ = σ (1 + ∆). Remember that σ representsthe softness of the agar medium.Regarding the chemotactic term, the coefficient χ represents the strengthof the signal. The simulations were performed using several values of this pa-rameter including the value χ = 0, that is, in the absence of chemotaxis.Thus, we also computed solutions to the original system (1) as well (in its non-dimensional form) also for the sake of comparison. For a complete summary ofthe free parameter values used during the numerical experiments, see table 1. We now present the numerical simulations of system (12) subject to initialconditions (15) and n ( x, y,
0) = n , as well as to boundary conditions (14). Wedivide the exposition of our numerical results in terms of the agar’s hardness.The latter is expressed by the parameter value σ . In our simulations, weconsidered the values σ = 1 and σ = 4, and the observed patterns correspondto regions E and D in the morphology diagram (see the enclosed region insidea rectangle in figure 1). For σ = 1 we took the time step ∆ t = 0 . σ = 4 we took it as ∆ t = 0 . χ , and each row is associated to a different time step.Such configuration allows the reader to compare the colony patterns for differentchemotactic intensities. In what follows we highlight that the main effect of thechemotaxis is the increment on the growth speed of the outer envelope of thebacterial colony. 9 .1.1. Semi-solid agar We begin with the case of semi-solid agar, for which σ = 1. This meansthat the concentration of agar is intermediate and the diffusion of the bacteria islimited. The observed patterns correspond to the DBM-like branching structurein region E of the morphology diagram. The considered values for the initialnutrient level were n = 0 .
71 and 1 . n = 0 .
71. On the left columnwe find the solutions to the system in the absence of chemotaxis, that is, when χ = 0. The pictures for this case are comparable to those obtained by Kawasaki et al. (see figure 3(b) in [1]). On the right column we find the case where thechemotactic sensitivity is set as χ = 5. The morphology shown in both casesis very similar (if we look carefully we can see that the number of branches isconsistent). There is, however, a significant increase in the growth velocity ofthe outer envelope of the branches when the chemotactic signal is switched on,as the reader may observe.The simulations for n = 1 .
07 are depicted in figure 6. Once again theleft column corresponds to the value χ = 0. The patterns are defined by fineradially oriented branches with round envelope. At the center of the colonythere is the highest concentration of bacteria. Actually, the peak concentrationis the remainder of the initial inoculation, because the bacteria at early timesdeplete the nutrient in site of the initial condition, but the cells in the innerpart of the colony become immotile. Meanwhile, the cells at the periphery ofthe colony move actively searching for higher concentrations of nutrient; this iswhat drives the formation of branches. Hence, bacteria at the front of the colonyaccount for the instability at the interface. This is the main feature of the non-linear cross diffusion model. In the right column of figure 6 we find the patternsin the presence of a chemotactic signal ( χ = 5 . In our next numerical experiment, we set the value σ = 4, causing thediffusion of the bacteria to be increased by lowering the concentration of theagar. At sufficient high levels of nutrient the colony patterns are restricted toregion D in the morphology diagram, whereas at lower levels of n patternsin region E emerge. In order to explore the morphology transitions betweenregions D and E, we first set n = 0 .
71 as the initial nutrient and considervarious values for χ , because the chemotactic signal gives rise to a outwarddrift to the cellular movement. The computations are depicted in figures 7and 8. In this experiment, the comparison is made through the two figures,for which the initial nutrient is fixed and the chemotactic sensitivity takes thevalues χ = 0 , . , . .
5. In figure 7, the left column shows the evolutionof the system without chemotaxis. At early times the bacteria evolve as a solid10olony (figure 7(a)), but for the consequent iterations they develop a patternof radially oriented fjords. On the right column of figure 7 we show the resultswith chemotactic strength χ = 2 .
5. As expected, the chemotactic term makesthe colony grow faster; compare, for example, figures 7(f) and 7(g). Observethat the pictures in the two columns of figure 7 (with chemotactic sensitivities χ = 0 and χ = 2 .
5) exhibit similar morphologies, being the envelope velocitythe only apparent difference between them.For higher values of χ , however, the growth speed is not only greatly in-creased, but there is also a change in the morphologies as well. We may appre-ciate, for instance, that for χ = 7 . t = 324 . χ the propagation speed of the front is higher than thediffusion of the nutrient preventing the fjords from forming. Therefore, the thinbranches start fusing (figure 8(g)) until the pattern shows no more branchesand becomes an homogeneous disk (figure 8(h)). These observations suggestthat, in the transition between regions E and D in the morphology diagram,the supressing effect of the chemotactic term on the onset of instability (whichcauses the formation of branches; cf. [30]) is more evident than in other regionsof the diagram. An increase of the chemotactic intensity prevents the formationof further branches and results in an homogeneization of the colony morphology,accompanied with the expected increase in the velocity of the overall colony.Finally, in figure 9 we present the results of our numerical simulations for thevalue of nutrient n = 1 .
07 and chemotaxis sensitivities χ = 0 and χ = 7 .
5. Inthis case no significant morphological changes are observed. In the left columnof figure 9 the pattern shown is an homogeneous disk with a peak of bacterialconcentration at the center. In the right column (when χ = 7 . It has been suggested that numerical simulations might depend on the gridsize [31]. To circumvent this possibility, we performed simulations with thefollowing grid sizes as well: N = 1024 , , σ = 1, n = 0 .
71 and χ = 5 can be observed in figure 2 foreach grid size. Notice that in the case of N = 1024 points per side, the solutionsexhibit a cross shaped pattern with a rhomboid envelop (figure 2(a)). Whenthe grid size is increased to N = 2048 there is indeed a significant change inmorphology, showing a denser branching pattern with a round envelope (figure2(b)). This pattern does not change significantly, however, when the size gridis increased to N = 4096 and N = 8192 (figures 2(c) and 2(d)). The readermay notice that figures 2(b), 2(c) and 2(d) show very similar patterns. We11onclude that a resolution of N = 2048 lies above the threshold needed for areliable display of the solutions of the system. Even though our computationswere performed with parallelization techniques on a GPU, the simulations with N = 8192 involved a large computational time (there are more than 64 milliongrid points). Therefore, we adopted N = 2048 as the resolution for all thesimulations of system (12) in this paper. (a) N = 1024 (b) N = 2048(c) N = 4096 (d) N = 8192Figure 2: Simulations of system (12), with parameter values σ = 1 , n = 0 .
71 and χ = 5,for different values of the grid size N .
4. Propagation velocity of the front
In this section we discuss the effects of the chemotactic term on the envelopepropagation speed by means of an asymptotic geometric front approximation.Moreover, we provide numerical estimations of the speed by simulating a onedimensional version of system (12). We compare these velocity calculations tothose without chemotaxis. 12 .1. Approximation of the speed
For simplicity, let us consider the chemotactic sensitivity function as a con-stant, i.e., χ = χ . Therefore, the system of equations takes the form n t = ∆ n − nb (16a) b t = ∇ · ( σnb ∇ b ) + nb − χ σ ∇ · (cid:0) nb ∇ n (cid:1) , (16b)In order to approximate the velocity of the smooth envelope enclosing thebranching patterns, we explore the following reduction proposed by Kawasaki etal. [1]. They made a further approximation in order to obtain a scalar equationfor the bacterial concentration b . They observed that in the absence of diffu-sion the total mass is conserved. In our setting, if we ignore both diffusion andchemotaxis then we can add the equations and integrate them in time to obtain n + b = C (17)where C is a constant, which will be taken as C = n . Thus, n = n − b and weobtain the logistic equation b t = n b (cid:18) − bn (cid:19) . It must be observed that in [1] the term σ nb is substituted by σ n b , andthe term nb is substituted by n b (1 − b/n ), in the approximation of the systemby a scalar equation for b . Instead, we are substituting (17) into (16b) to obtaina scalar equation for b . This not only represents a more accurate approximationbut it is also consistent with the numerical computations which, at first order,show that the conservation-like equation (17) is valid away from the front.Substitution of (17) into (16) yields the scalar equation b t = ∇ · ( σ ( n − b ) b ∇ b ) + ( n − b ) b − σ χ ∇ · (( n − b ) b ∇ ( n − b )) , (18)which is an approximated scalar equation for b with a Fisher-type reactionterm. Here we have approximated σ (1 + ∆) ≈ σ , because in this part of thestudy we are not interested in considering the random motion of bacteria. Arearrangement of the terms in equation (18) allows us to write the two divergenceterms together. This yields b t = σ ∇ · ( ˜ D ( b, χ ) ∇ b ) + g ( b ) , (19)where ˜ D ( b, χ ) = n b (cid:18) − bn (cid:19) (1 + χ b ) , (20)may be interpreted as an effective non-linear diffusion coefficient, which involvesboth the contributions of diffusion and chemotaxis. The reaction term is ofFisher-KPP type, g ( b ) = n b (cid:18) − bn (cid:19) . (21)13e describe the motion of the front in local curvilinear coordinates withcomponents ζ ( x, y, t ) and τ ( x, y, t ), the normal and tangent unit vectors to thefront, respectively. According to custom, ζ is the normal component pointinginside the colony front, so that − ζ t is the outer normal speed. We assume thatthe dependence of b with respect to the tangent component is negligible and weapproximate b ( x, y, t ) ≈ ¯ b ( ζ ( x, y, t )) . (22)After substituting ¯ b into equation (19), and dropping the bar for notation con-venience, we get the following ordinary differential equation for b − cb (cid:48) = σ ( ˜ D ( b ) b (cid:48) ) (cid:48) + b ( n − b ) . (23)where c = − ζ t + σ ˜ D ( b )∆ ζ (24)and − ζ t is the velocity of the envelope front that we want to estimate.The approximated front equation (23) in the normal direction resembles aone-dimensional equation of form (A.6) (see Appendix A). By the results ofMalaguti and Marcelli [32] the one-dimensional velocity for each value of χ ≥ χ , n , and the functions ˜ D and g . We note that weare interested in comparing the velocity between the cases χ = 0 and χ > c ( χ ) bounded by c ( χ ) ≥ c ∗ ( χ ), where the threshold velocity c ∗ satisfies the bound:0 < c ∗ ( χ ) ≤ ¯ c ( n , χ ) := 2 (cid:115) sup b ∈ (0 ,n ] σ ˜ D ( b, χ ) g ( b ) b (25)with ˜ D and g given by (20) and (21). Upon substitution, we are able to computethe bound ¯ c ( n , χ ) = 2 n √ σ (cid:115) max b ∈ [0 ,n ] b (cid:18) − bn (cid:19) (1 + χ b ) . (26)The function ˜ ψ ( b ) = b (1 − b/n ) (1 + χ b ) is non-negative for 0 ≤ b ≤ n and each χ ≥
0; moreover, it has a global maximum ˜ ψ max = ˜ ψ ( b ∗ ), where b ∗ = 12 (cid:18) n − χ (cid:19) + 12 (cid:115)(cid:18) n − χ (cid:19) + n χ , (27)and with 0 < b ∗ < n for all values of χ > n > c ( n , χ ) = 2 n √ σ (cid:112) ψ ( b ∗ ) . (28)In the absence of chemotaxis ( χ = 0) the function ˜ ψ ( b ) reduces to ψ ( b ) = b (1 − b/n ) and has a global maximum ψ ( n /
3) = 4 n /
27. This means that¯ c ( n ,
0) = 2 n √ σ (cid:114) n
27 = 43 √ n / σ / . (29)14quations (28) and (29) are the speed thresholds for the cases with and withoutchemotaxis, respectively. These equations predict the speed of a sharp front, andwill be used when comparing the numerical estimates of the velocity (see section4.2 below). In our setting, we are considering that the growth of the colony is inthe radial direction (see equation (22)), so we only consider the normal velocity s = − ζ t . Furthermore, we take into account the cases where the colony exhibitsan outer growth, this means that the local curvature is positive there κ = − ∆ ζ > . (30)Returning to equations (23) and (24) and comparing the one dimensionalvelocity of a front for equation (A.6) to expression (24), we obtain s = − ζ t ≥ c ∗ ( χ ) + σ ˜ D ( b ) κ > . (31)It follows from equation (30) and Proposition 1, in Appendix A, that s = − ζ t ≥ c ∗ ( χ ) + σ ˜ D ( b ) κ ≥ c ∗ (0) + σ ˜ D ( b ) κ ≥ c ∗ (0) . (32)Last inequality implies that the normal velocity s is greater than the ve-locity for the sharp front when there is no chemotaxis ( χ = 0), and that itis increased/decreased by a term proportional to the local curvature and thechemotactic strength χ >
0. In the outer front where diffusion is degenerate,this last curvature term is negligible (as ˜ D = 0 in a vicinity of the envelopefront) so that we can approximate the speed by the one-dimensional sharp frontvelocity, c ∗ ( χ ) (cid:46) s , regardless of the curvature sign. As we have stated in (22), the tangential direction is negligible in describingthe front movement. Then, we shall estimate the velocity of the envelope frontsolving numerically a one-dimensional version of system (16), that is of thefollowing form, ∂n∂t = (cid:18) ∂n∂x (cid:19) x − nb (33a) ∂b∂t = (cid:18) σ nb ∂b∂x − χ σ nb ∂n∂x (cid:19) x + nb. (33b)The simulations of system (33) were performed on the spacial domain [0 , n ( x,
0) = n , b ( x,
0) = b M e − x / . , (34)where b M = 0 .
71. Again, as in the two-dimensional simulations, n measures theinitial level of nutrient, χ measures the chemotaxis intensity and σ measuresthe agar hardness. We made several simulations of system (33) for variousvalues of n = { . , . , . , , , } , χ = { , , . , , } . In all simulations15e fixed σ = 1 . Figure 3 depicts the numerical solutions of the bacterial densityof b at different time steps and different values of n and χ . It should be notedthat solutions of b initially do not exhibit a traveling wave profile, nevertheless,eventually the solution resembles a traveling wave. Therefore, the traveling waveestimations were computed at later times by selecting a point in the vicinity of0.5. These are the points shown in figure 3. In figure 4 we present a comparisonof the propagation velocity of the front between the numerical estimations andthe velocity threshold given by equations (28) for χ > χ = 0, asfunctions of the initial nutrient n . For the case without chemotaxis, we foundsharper results (see figure 4(a)) than those in the original work by Kawasaki et al. , because we are considering a more accurate velocity threshold, as wasexplained in section 4.1. As expected from equation (32), when a chemotacticsignal is present the velocity of the front is greatly increased. We can see thatthe numerical velocity fits the velocity threshold. Therefore, figure 4 highlightstwo facts: first, numerical estimations of the propagation speed of the one-dimensional system are in good accordance with the theoretical speed thresholdsthat are defined for the scalar equation for b (see equation (19)); and second,the conservation-like equation (17) is a good approximation of the solutions ofsystem (16).
5. Conclusions
In this work we have studied the effects that nutient chemotaxis inducesinto bacterial population patterns modeled by a reaction-diffusion system orig-inally proposed by Kawasaki et al. [1]. Motivated by a discussion provided byGolding et al. [6], we introduce into system (1) a chemotactic term compatiblewith the non-linear cross diffusion. After a suitable rescaling of variables, wereduce the number of free parameters of system (9), and we only considered σ , n and χ . We solved numerically system (12), for various values of the freeparameters, including the case without chemotaxis ( χ = 0). This means thatwe reproduced the numerical results obtained by Kawasaki et al. in [1]. Ournumerical simulations of system (12) show that the chemotactic term providesan outward drift to the bacterial movement and therefore enhances the growthvelocity of the bacterial patterns.In order to approximate the velocity of the envelope front, under suitableassumptions, we derived a scalar equation of the bacterial density, see equation(19). Then, applying a result from Malaguti and Marcelli [32], see AppendixA, we proved that the normal velocity s is greater than the velocity of the sharpfront when there is no chemotaxis ( χ = 0). Afterwards, by simulating a one-dimensional version of system (16), we calculated numerical estimations of thepropagation velocity of the front. We compared the theoretical speed predic-tions, defined in equations (28) and (29), versus the estimations and we foundthat the propagation velocity of the envelope front is greatly increased when achemotactic signal is present (see figure 4). For the case without chemotaxis wefound results comparable with those found by Kawasaki et al. (see figure 4(a)),although we must stress the fact that we are using a more accurate velocity16 a) n = 0 . χ = 10 (b) n = 2, χ = 2 . n = 3, χ = 1 (d) n = 4, χ = 1Figure 3: Numerical solution of the bacterial density b from the system (33) at different timesand different values of n and χ . threshold than the one used in [1]. Furthermore, as figure 4 depicts, the nu-merical estimations of the propagation speed of the one-dimensional system arein good accordance with the theoretical speed thresholds that are defined forthe scalar equation for b (see equation (19)), implying that the conservation-likeequation (17) is a good approximation of the solutions of system (16).Finally, we discuss the relation of our results to the analysis of Arouh andLevine [30], who studied the Kessler-Levine reaction-diffusion equations (cf.[33]) and explored the addition of a nutrient chemotactic term. They foundthat the chemotaxis suppresses the onset of instability causing the formationof branching patterns in the colony envelope as it propagates outward. Thissuppression of instability is reflected in a decrease of the “growth rate” of theperturbation of the envelope when the chemotaxis is switched on. The au-thors conjecture that this behavior can be extrapolated to other (more realistic)chemotactic models like the one introduced in this paper. It is to be pointedout, however, that our findings do not contradict the analysis of Arouh andLevine inasmuch as our numerical simulations and analytic computations ex-hibit an increase in the normal velocity of the envelope outward due to the17 a) χ = 0 (b) χ = 1 . χ = 2 . χ = 5 . χ = 10 . χ . The asterisks depict the numerical estimations of the velocity. The solid lineis the velocity threshold given by equation (29) in the absence of chemotaxis (case 4(a)), andby equation (28) when chemotaxis is switched on (cases 4(b), 4(c), 4(d) and 4(e)), both asfunctions of the initial nutrient concentration n . Here σ = 1. chemotaxis, whereas the growth rate that Arouh and Levine refer to is the tem-poral eigenvalue associated to an eigenfunction (perturbation) of the linearized18perator around the envelope front. Actually, a simple calculation using energyestimates (not included here) shows that in the case of our model, the spectralbounds for the eigenvalues of the linearized operator around the envelope frontindeed decrease as functions of the chemotactic sensitivity. Our results and theobservations of Arouh and Levine combined suggest that, although the speed ofthe envelope front increases with the chemotaxis (as it is shown in this paper),the growth rate of linear perturbations of the envelope is pushed to the stable(negative) region of the spectrum, delaying the formation of fingering unstablepatterns (as it is shown in [30]). Appendix A. Bounds for the velocity
In this section we summarize the theoretical results of Malaguti and Marcelli[32] that allow us to justify why the speed of the envelope front is greater whenthere is a chemotactic signal. Theorem 9 in [32] states that for one-dimensionalreaction-diffusion equations of the form u t = ( D ( u ) u x ) x + g ( u ) (A.1)with D ∈ C ([0 , g ∈ C ([0 , u = 0 and u = 1, namely, D (0) = D (1) = 0 , D ( u ) > < u < , (A.2)and the rate of growth of the population is of Fisher-KPP type, g (0) = g (1) = 0 , g ( u ) > < u < , (A.3)then there exist a traveling wave solution (unique up to translation) to equation(A.1), provided that the velocity c satisfies 0 < c ∗ ≤ c , where c ∗ is a thresholdvelocity satisfying the bound0 < c ∗ ≤ (cid:115) sup s ∈ (0 , D ( s ) g ( s ) s . (A.4)The case c = c ∗ corresponds to a sharp front. They prove the theorem byshowing that the existence of such a traveling wave is equivalent to the existenceof a solution z = z ( u ) of the boundary-value problem dzdu = − c − D ( u ) g ( u ) z ,z (0 + ) = z (1 − ) = 0 , (A.5)for the the same c , satisfying z ( u ) < u ∈ (0 , b t = ( ˜ D ( b, χ ) b x ) x + g ( b ) , (A.6)19n [0 , n ] and where ˜ D ( b, χ ) and g ( b ) are defined in equations (20) and (21),respectively. The existence of a traveling front solution to the equation with χ > dzdb = − c − ˜ D ( b, χ ) g ( b ) z ,z (0 + ) = z ( n − ) = 0 . (A.7)Likewise the traveling wave without chemotaxis is linked to the solutions to theproblem dzdb = − c − ˜ D ( b, g ( b ) z ,z (0 + ) = z ( n − ) = 0 . (A.8)In order to compare the speed of the envelope front in both cases, we considerthe following sets A χ = { c > } , (A.9) A = { c > } , (A.10)that define the range of values of c for which there exists a traveling wavesolution for equation (19) (see Theorem 9 in [32]). It should be noted that theinfimum of the set, in either case, corresponds to the threshold speed c ∗ . Thefollowing proposition compare the two velocity thresholds (with and withoutchemotaxis). Proposition 1.
For any fixed value of n > and for all values of χ ≥ therehold,(i) ¯ c ( χ ) ≥ ¯ c (0) , and(ii) c ∗ ( χ ) ≥ c ∗ (0) .Proof. The proof of (i) is straightforward: for all values of χ ≥ b ∈ (0 , n ) there holds ψ ( b, χ ) = b (1 − b/n ) (1 + χ b ) ≥ ψ ( b,
0) = b (1 − b/n ) . Therefore max b ∈ [0 ,n ] ψ ( b, χ ) ≥ max b ∈ [0 ,n ] ψ ( b, c > c ∈ A χ forsome χ >
0. This means that the problem (A.7) has a non-positive solution ζ ( b ). Therefore, since ˜ D ( b, χ ) > ˜ D ( b,
0) for all b ∈ (0 , n ) inasmuch as χ > dζdb = − c − ˜ D ( b, χ ) ζ ( b ) > − c − ˜ D ( b, ζ ( b ) . By Lemma 8 in [32], there exists a negative solution z = z ( b ) to problem (A.8),with the same constant c . This shows that c ∈ A . Therefore, A χ ⊆ A forall χ ≥
0, and consequently c ∗ ( χ ) = inf A χ ≥ c ∗ (0) = inf A . This yields(ii). 20 cknowledgements The work of J. F. Leyva was partially supported by CONACyT (Mexico)through a scholarship for doctoral studies, grant no. 220865.
References [1] K. Kawasaki, A. Mochizuki, M. Matsushita, T. Umeda, N. Shigesada, Mod-eling spatio-temporal patterns generated by
Bacillus subtilis , J. of Theor.Biol. 188 (2) (1997) 177 – 185.[2] T. Matsuyama, M. Matsushita, Fractal morphogenesis by a bacterial cellpopulation, Crit. Rev. Microbiol. 19 (2) (1993) 117–135.[3] E. Ben-Jacob, O. Schochet, A. Tenenbaum, I. Cohen, A. Czir´ok, T. Vic-sek, Generic modelling of cooperative growth patterns in bacterial colonies,Nature 368 (1994) 46–49.[4] M. Ohgiwari, M. Matsushita, T. Matsuyama, Morphological changes ingrowth phenomena of bacterial colony patterns, J. Phys. Soc. Jpn. 61 (3)(1992) 816–822.[5] J.-I. Wakita, K. Komatsu, A. Nakahara, T. Matsuyama, M. Matsushita,Experimental investigation on the validity of population dynamics ap-proach to bacterial colony formation, J. Phys. Soc. Jpn. 63 (3) (1994)1205–1211.[6] I. Golding, Y. Kozlovsky, I. Cohen, E. Ben-Jacob, Studies of bacterialbranching growth using reaction-diffusion models for colonial development,Phys. A 260 (3–4) (1998) 510–554.[7] M. Wu, X.-W. Zou, Z.-Z. Jin, Morphological diagram of bacterial colonypattern simulated using a revised Mimura-Sakaguchi-Matsushita model, J.Phys. Soc. Jpn. 74 (8) (2005) 2379–2380.[8] J. D. Murray, Mathematical biology II. Spatial models and biomedical ap-plications, 3rd Edition, Vol. 18 of Interdisciplinary Applied Mathematics,Springer-Verlag, New York, 2003.[9] J. A. Sherratt, On the form of smooth-front travelling waves in a reaction-diffusion equation with degenerate nonlinear diffusion, Math. Model. Nat.Phenom. 5 (5) (2010) 64–79.[10] R. A. Satnoianu, P. K. Maini, F. S. Garduno, J. P. Armitage, Travellingwaves in a nonlinear degenerate diffusion model for bacterial pattern for-mation, Discrete Contin. Dyn. Syst. Ser. B 1 (3) (2001) 339–362.[11] E. F. Keller, L. A. Segel, Model for chemotaxis, J. Theor. Biol. 30 (2)(1971) 225 – 234. 2112] E. F. Keller, L. A. Segel, Traveling bands of chemotactic bacteria: A the-oretical analysis, J. Theor. Biol. 30 (2) (1971) 235 – 248.[13] W. Alt, Biased random walk models for chemotaxis and related diffusionapproximations, J. Math. Biol. 9 (2) (1980) 147–177.[14] C. S. Patlak, Random walk with persistence and external bias, Bull. Math.Biophys. 15 (1953) 311–338.[15] T. Hillen, K. J. Painter, A user’s guide to PDE models for chemotaxis, J.Math. Biol. 58 (1-2) (2009) 183–217.[16] I. Cohen, A. Czir´ok, E. Ben-Jacob, Chemotactic-based adaptive self orga-nization during colonial development, Phys. A 233 (3–4) (1996) 678–698.[17] E. Ben-Jacob, From snowflake formation to the growth of bacterial coloniesII: Cooperative formation of complex colonial patterns, Contemp. Phys.38 (3) (1997) 205–241.[18] G. H. Wadhams, J. P. Armitage, Making sense of it all: bacterial chemo-taxis, Nature Rev. Molecular Cell Biol. 5 (12) (2004) 1024–1037.[19] M. Eisenbach, Bacterial chemotaxis, eLS. John Wiley and Sons, Chichester(2011).[20] R. I. Lapidus, R. Schiller, Model for the chemotactic response of a bacterialpopulation, Biophys. J. 16 (7) (1976) 779–789.[21] R. Ford, D. Lauffenburger, Measurement of bacterial random motility andchemotaxis coefficients: II. Application of single-cell-based mathematicalmodel, Biotechn. Bioeng. 37 (7) (1991) 661–672.[22] L. A. Segel, Incorporation of receptor kinetics into a model for bacterialchemotaxis, J. Theoret. Biol. 57 (1) (1976) 23–42.[23] L. A. Segel, A theoretical study of receptor mechanisms in bacterial chemo-taxis, SIAM J. Appl. Math. 32 (3) (1977) 653–665.[24] R. Mesibov, G. W. Ordal, J. Adler, The range of attractant concentrationsfor bacterial chemotaxis and the threshold and size of response over thisrange, J. Gen. Physiol. 62 (2) (1973) 203–223.[25] H. G. Othmer, A. Stevens, Aggregation, blowup, and collapse: the ABCsof taxis in reinforced random walks, SIAM J. Appl. Math. 57 (4) (1997)1044–1081.[26] D. J. Goldman, G. W. Ordal, Sensory adaptation and deadaptation by bacillus subtilis , J. of Bacteriology 147 (1) (1981) 267–70.[27] S. Kitsunezaki, Interface dynamics for bacterial colony formation, J. Phys.Soc. Japan 66 (5) (1997) 1544–1550.2228] D. B. Kirk, W.-M. W. Hwu, Programming Massively Parallel Processors:A Hands-on Approach, 1st Edition, Morgan Kaufmann Publishers Inc., SanFrancisco, CA, USA, 2010.[29] NVIDIA, NVIDIA CUDA C Programming Guide 4.0, 2011.[30] S. Arouh, H. Levine, Nutrient chemotaxis suppression of a diffusive insta-bility in bacterial colony dynamics, Phys. Rev. E 62 (1) (2000) 1444–1447.[31] M. Mimura, H. Sakaguchi, M. Matsushita, Reaction diffusion modelling ofbacterial colony patterns, Phys. A 282 (12) (2000) 283–303.[32] L. Malaguti, C. Marcelli, Sharp profiles in degenerate and doubly degener-ate Fisher-KPP equations, J. Differential Equations 195 (2) (2003) 471–496.[33] D. A. Kessler, H. Levine, Fluctuation-induced diffusive instabilities, Nature394 (1998) 556–558. 23 a) χ = 0, t = 694 . χ = 5, t = 694 . χ = 0, t = 1389 . χ = 5, t = 1389 . χ = 0, t = 2083 . χ = 5, t = 2083 . χ = 0, t = 2824 . χ = 5, t = 2824 . σ = 1 . , n = 0 . χ = 0 (no chemotaxis - left), and χ = 5 . t . a) χ = 0, t = 115 . χ = 5, t = 115 . χ = 0, t = 231 . χ = 5, t = 231 . χ = 0, t = 370 . χ = 5, t = 370 . χ = 0, t = 509 . χ = 5, t = 509 . σ = 1 . , n = 1 . χ = 0 (no chemotaxis - left), and χ = 5 . t . a) χ = 0, t = 92 . χ = 2 . t = 92 . χ = 0, t = 162 . χ = 2 . t = 162 . χ = 0, t = 231 . χ = 2 . t = 231 . χ = 0, t = 324 . χ = 2 . t = 324 . σ = 4 . , n = 0 . χ = 0 (left) and χ = 2 . t ..
Bacillus subtilis , J. of Theor.Biol. 188 (2) (1997) 177 – 185.[2] T. Matsuyama, M. Matsushita, Fractal morphogenesis by a bacterial cellpopulation, Crit. Rev. Microbiol. 19 (2) (1993) 117–135.[3] E. Ben-Jacob, O. Schochet, A. Tenenbaum, I. Cohen, A. Czir´ok, T. Vic-sek, Generic modelling of cooperative growth patterns in bacterial colonies,Nature 368 (1994) 46–49.[4] M. Ohgiwari, M. Matsushita, T. Matsuyama, Morphological changes ingrowth phenomena of bacterial colony patterns, J. Phys. Soc. Jpn. 61 (3)(1992) 816–822.[5] J.-I. Wakita, K. Komatsu, A. Nakahara, T. Matsuyama, M. Matsushita,Experimental investigation on the validity of population dynamics ap-proach to bacterial colony formation, J. Phys. Soc. Jpn. 63 (3) (1994)1205–1211.[6] I. Golding, Y. Kozlovsky, I. Cohen, E. Ben-Jacob, Studies of bacterialbranching growth using reaction-diffusion models for colonial development,Phys. A 260 (3–4) (1998) 510–554.[7] M. Wu, X.-W. Zou, Z.-Z. Jin, Morphological diagram of bacterial colonypattern simulated using a revised Mimura-Sakaguchi-Matsushita model, J.Phys. Soc. Jpn. 74 (8) (2005) 2379–2380.[8] J. D. Murray, Mathematical biology II. Spatial models and biomedical ap-plications, 3rd Edition, Vol. 18 of Interdisciplinary Applied Mathematics,Springer-Verlag, New York, 2003.[9] J. A. Sherratt, On the form of smooth-front travelling waves in a reaction-diffusion equation with degenerate nonlinear diffusion, Math. Model. Nat.Phenom. 5 (5) (2010) 64–79.[10] R. A. Satnoianu, P. K. Maini, F. S. Garduno, J. P. Armitage, Travellingwaves in a nonlinear degenerate diffusion model for bacterial pattern for-mation, Discrete Contin. Dyn. Syst. Ser. B 1 (3) (2001) 339–362.[11] E. F. Keller, L. A. Segel, Model for chemotaxis, J. Theor. Biol. 30 (2)(1971) 225 – 234. 2112] E. F. Keller, L. A. Segel, Traveling bands of chemotactic bacteria: A the-oretical analysis, J. Theor. Biol. 30 (2) (1971) 235 – 248.[13] W. Alt, Biased random walk models for chemotaxis and related diffusionapproximations, J. Math. Biol. 9 (2) (1980) 147–177.[14] C. S. Patlak, Random walk with persistence and external bias, Bull. Math.Biophys. 15 (1953) 311–338.[15] T. Hillen, K. J. Painter, A user’s guide to PDE models for chemotaxis, J.Math. Biol. 58 (1-2) (2009) 183–217.[16] I. Cohen, A. Czir´ok, E. Ben-Jacob, Chemotactic-based adaptive self orga-nization during colonial development, Phys. A 233 (3–4) (1996) 678–698.[17] E. Ben-Jacob, From snowflake formation to the growth of bacterial coloniesII: Cooperative formation of complex colonial patterns, Contemp. Phys.38 (3) (1997) 205–241.[18] G. H. Wadhams, J. P. Armitage, Making sense of it all: bacterial chemo-taxis, Nature Rev. Molecular Cell Biol. 5 (12) (2004) 1024–1037.[19] M. Eisenbach, Bacterial chemotaxis, eLS. John Wiley and Sons, Chichester(2011).[20] R. I. Lapidus, R. Schiller, Model for the chemotactic response of a bacterialpopulation, Biophys. J. 16 (7) (1976) 779–789.[21] R. Ford, D. Lauffenburger, Measurement of bacterial random motility andchemotaxis coefficients: II. Application of single-cell-based mathematicalmodel, Biotechn. Bioeng. 37 (7) (1991) 661–672.[22] L. A. Segel, Incorporation of receptor kinetics into a model for bacterialchemotaxis, J. Theoret. Biol. 57 (1) (1976) 23–42.[23] L. A. Segel, A theoretical study of receptor mechanisms in bacterial chemo-taxis, SIAM J. Appl. Math. 32 (3) (1977) 653–665.[24] R. Mesibov, G. W. Ordal, J. Adler, The range of attractant concentrationsfor bacterial chemotaxis and the threshold and size of response over thisrange, J. Gen. Physiol. 62 (2) (1973) 203–223.[25] H. G. Othmer, A. Stevens, Aggregation, blowup, and collapse: the ABCsof taxis in reinforced random walks, SIAM J. Appl. Math. 57 (4) (1997)1044–1081.[26] D. J. Goldman, G. W. Ordal, Sensory adaptation and deadaptation by bacillus subtilis , J. of Bacteriology 147 (1) (1981) 267–70.[27] S. Kitsunezaki, Interface dynamics for bacterial colony formation, J. Phys.Soc. Japan 66 (5) (1997) 1544–1550.2228] D. B. Kirk, W.-M. W. Hwu, Programming Massively Parallel Processors:A Hands-on Approach, 1st Edition, Morgan Kaufmann Publishers Inc., SanFrancisco, CA, USA, 2010.[29] NVIDIA, NVIDIA CUDA C Programming Guide 4.0, 2011.[30] S. Arouh, H. Levine, Nutrient chemotaxis suppression of a diffusive insta-bility in bacterial colony dynamics, Phys. Rev. E 62 (1) (2000) 1444–1447.[31] M. Mimura, H. Sakaguchi, M. Matsushita, Reaction diffusion modelling ofbacterial colony patterns, Phys. A 282 (12) (2000) 283–303.[32] L. Malaguti, C. Marcelli, Sharp profiles in degenerate and doubly degener-ate Fisher-KPP equations, J. Differential Equations 195 (2) (2003) 471–496.[33] D. A. Kessler, H. Levine, Fluctuation-induced diffusive instabilities, Nature394 (1998) 556–558. 23 a) χ = 0, t = 694 . χ = 5, t = 694 . χ = 0, t = 1389 . χ = 5, t = 1389 . χ = 0, t = 2083 . χ = 5, t = 2083 . χ = 0, t = 2824 . χ = 5, t = 2824 . σ = 1 . , n = 0 . χ = 0 (no chemotaxis - left), and χ = 5 . t . a) χ = 0, t = 115 . χ = 5, t = 115 . χ = 0, t = 231 . χ = 5, t = 231 . χ = 0, t = 370 . χ = 5, t = 370 . χ = 0, t = 509 . χ = 5, t = 509 . σ = 1 . , n = 1 . χ = 0 (no chemotaxis - left), and χ = 5 . t . a) χ = 0, t = 92 . χ = 2 . t = 92 . χ = 0, t = 162 . χ = 2 . t = 162 . χ = 0, t = 231 . χ = 2 . t = 231 . χ = 0, t = 324 . χ = 2 . t = 324 . σ = 4 . , n = 0 . χ = 0 (left) and χ = 2 . t .. a) χ = 5, t = 92 . χ = 7 . t = 92 . χ = 5, t = 162 . χ = 7 . t = 162 . χ = 5, t = 231 . χ = 7 . t = 231 . χ = 5, t = 324 . χ = 7 . t = 324 . σ = 4 . , n = 0 . χ = 5 . χ = 7 . t . a) χ = 0, t = 33 . χ = 7 . t = 33 . χ = 0, t = 66 . χ = 7 . t = 66 . χ = 0, t = 88 . χ = 7 . t = 88 . χ = 0, t = 110 . χ = 7 . t = 110 . σ = 4 . , n = 1 . χ = 0 (left) and χ = 7 . t ..