Computational Studies of Bacterial Colony Model

Microbiological experiments show that the colonies of the bacterium bacillus subtilis placed on a dish filled with an agar medium and nutrient form varied spatial patterns while the individual cells grow, reproduce and migrate on the dish in clumps. In this paper, we discuss a system of reaction-diffusion equations that can be used with a view to modelling this phenomenon and we solve it numerically by means of the method of lines. For the spatial discretization, we use the finite difference method and Galerkin finite element method. We present how the spatial patterns obtained depend on the spatial discretization employed and we measure the experimental order of convergence of the numerical schemes used. Further, we present the numerical results obtained by solving the model in a cubic domain.


Introduction
According to microbiological experiments, the colonies of the bacterium bacillus subtilis placed on a dish filled with an agar medium and nutrient form varied-and surprisingly regular-patterns while the individual cells grow, reproduce and migrate on the dish in clumps (see, e.g., [1,2]).As bacterium is an extremely primitive organism, the question arises whether the motion of the cells is governed by some simple underlying principles that can be described by mathematically.
The general approach to the pattern formation in biology is based on continuous and discrete dynamical models or neural models as described in [3].In this case, the simple character of the motion of the bacteria resembles the Brownian motion and indicates the analogy with the diffusion limited aggregation (see [4]).Detailed analysis of the results of the experiments coupled the diffusion effects with the reactive behaviour of the species and led to the derivation of the model based on the reaction-diffusion system (see [5,6]) which was studied in [2,7].In this article, we discuss it with respect to the accuracy and reliability of the numerical solution.
The model is based on the assumption that the cells can be divided into two classes: the active and the inactive.The active cells, with concentration at the point x and time t, grow, move, feed on the nutrient and reproduce; whereas the inactive ones, with concentration  (1) for x   and , where d represents a diffusion coefficient, and the function is of the form with positive constants 0 1 .The system (1)-( 3) is subject to the initial conditions 2 , , a a a  from the Sobolev space in the sense of distributions on a time interval , where stands for the inner product , where ; for all x   in the sense of distributions on I ; and where The existence and uniqueness of the weak solution are proved in [8].The authors also proved that for a positive constant , and that there is a function (a final colony pattern) such that terns in 2    (the darker points correspond to the higher concentration of the cells) that have been found by the author by solving system (1)-(3).

Numerical Schemes
The model under consideration was solved numerically by means of the method of lines; for the spatial discretization, the finite difference method and Galerkin finite element method (the latter only in the case 2    ) were used.For the sake of simplicity, we use only square and cubic domains  , specifically

Finite Difference Method
In the two-dimensional case, the domain is covered by a regular orthogonal mesh of nodes are defined similarly).The spatial derivatives in (1) and ( 2) are approximated by the second order central difference, and the derivatives in the boundary conditions are approximated by the first order backward difference.Therefore, we obtain the following system of ODE's: , , , where in the first two equations, ; and in the third equation, Further, on the boundary we get , ; This system is supplemented with the initial conditions , In the case of the cubic domain , the spatial discretization is derived analogically.These two schemes will be called MOL.

Galerkin Finite Element Method
In this method, both the linear (triangular) and the bilinear (square) Lagrange elements (see [9]) are used.The basis   W  consists of the functions which are linear on each triangle (in the first case) or bilinear on each square (in the second case) and take the value 1 at one node of the spatial mesh and vanish at the other nodes.In the first case, two types of triangulations are used; they are depicted in Figure 2. In the case of the square finite elements, the same spatial mesh as in the method of lines is utilized.Further, the integrand in 5) is approximated by its interpolant (see [9]).Finally, the numerical scheme is simplified employing the method of lumped masses (see [10,11]).Hence, writing and , where N h is the number of the nodes of the spatial mesh, Equations ( 5)-( 7) are transformed as follows: , where h I stands for the interpolant.
This system is supplemented with the initial conditions The schemes corresponding to the cases depicted in Figure 2 will be denoted FEMHS, FEMRHS and FEMSS, respectively.
The convergence of the numerical schemes derived can be proved as in [12][13][14].The systems of ODE's derived were solved by means of the Runge-Kutta-Merson method (see [12,15]) with the adaptive time step control.

Quantitative Studies
In the two-dimensional case, the experimental order of convergence (EOC), see [16], was measured; only the error of the functions u and v was considered.As the model studied does not possess an analytical solution, Equations ( 1) and ( 2) were transformed into the form where the additional terms
The modified system was solved analogically as the original one (the functions f and g were approximated in the same way as the functions u and v in ( 9)).The EOC coefficient  is defined by where , i h represents the numerical solution corresponding to the parameter i which denotes the constant distance between two nodes on the bottom side (according to Figure 2) of the spatial mesh used (there are N nodes on this side), and  denotes a suitable norm.The following form of Equation ( 14) was used: where stands for a final time, and T  is a time step between the compared time levels of the solution.Two norms where denote the nodal values of the projections of and onto a regular orthogonal mesh with a spatial step .

Computational Studies of Pattern Formation
This section provides the comparison of final spatial patterns obtained by the four schemes under consideration in a square domain.All the results were computed for The results presented (Figures 3-11) have the form of coloured pictures with a final shape of the function w, computed by one of the schemes, in shades of grey; each of the coloured lines in these pictures traces the spatial shape of the function w computed by one of the schemes for one value of N. Further, in Tables 5-7 there are differences between these shapes expressed by means of the -and -norm (plus one another table for the combination 0 , ).The value of the 2norm is divided by the area of .The time evolution of the patterns is shown in [17].
Note that, in general, it cannot be stated that one scheme is more or less accurate than the other schemes because, typically, each of them produces another shape.In author's opinion, the most significant difference between the schemes consists in the symmetry of the results obtained by the scheme FEMHS.

Conclusions
It follows from the numerical results that the shape of the numerical solution strongly depends on the spatial discretization employed.The pattern type, nevertheless, seems to be stable, if the spatial mesh is not too coarse.In author's opinion, therefore, the numerical methods under consideration can be used for the investigation of the pattern growth described by the system of the reaction-diffusion equations discussed.
As for the measurement of the experimental order of convergence, it is remarkable that the method of lumped masses does not seem to decrease it.Finally, the threedimensional model exhibits very complex behaviour as well.

Acknowledgements
Partial support of the projects "Applied Mathematics in Physics and Technical Sciences" number MSM684077 0010 of the Ministry of Education, Youth and Sports of the Czech Republic and "Advanced Supercomputing Methods for Implementation of Mathematical Models" number SGS11/161/OHK4/3T/14 is acknowledged.
do anything.Finally, let be the concentration of the nutrient.
Then the time and space evolution of the concentrations within a bounded spatial domain , with  or 3 n  , is governed by the system

Figure 1 Figure 1 .
Figure1shows the observed types of the spatial pat-

Figure 2 .
Figure 2. Meshes used for the Galerkin FEM.
chosen constant.Therefore, the pattern types vary depending on the values of the parameters and .v d 0

Figure 10 .
Figure 10.Comparison of the patterns computed by the scheme MOL for v 0 = 0.065, d = 0.1 and different values of N.

Figure 11 .Table 5 .
Figure 11.Comparison of the patterns computed by the scheme FEMHS for v 0 = 0.1, d = 0.05 and different values of N.
This section provides examples of the spatial patterns that arise if we consider the cubic domain with 200 L  or .Similar types of the initial conditions and values of the parameters as in the two-dimensional case 300 L  were used, i.e., third spatial coordinate).The values of 0 and d represented variable parameters of the model.The results are presented in Figures 12-25.They have the form of coloured isolines of the functions u, v or .The time evolution is shown only for one combination of parameters; in the other cases, the solution evolves analogically.

Table 1 . Results of the EOC measurement for the scheme MOL.
refThe results are summarized in Tables1-4.The EOC coefficients in the i-th row are computed from the data in the i-th and -th row.The schemes FEMRHS, FEMHS and FEMSS exhibit the second order of convergence in h, and MOL exhibits the first order.Further,

Table 4 . Results of the EOC measurement for the scheme FEMRHS.
FEMSS appears to be the most accurate of all.