ABSTRACT

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.

1. 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, do not 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, is governed by the system

(1)

(2)

(3)

for and, where d represents a diffusion coefficient, and the function is of the form

with positive constants. The system (1)-(3) is subject to the initial conditions

(4)

for where are functions (represents the closure of), and is a positive constant; and the boundary conditions

where denotes the boundary of, and stands for the derivative of in the direction of the unit outward normal.

This problem can be written in the following weak formulation:

(5)

(6)

for all from the Sobolev space in the sense of distributions on a time interval, where stands for the inner product

, where;

(7)

for all in the sense of distributions on; and

(8)

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

Figure 1 shows the observed types of the spatial pat-

terns in (the darker points correspond to the higher concentration of the cells) that have been found by the author by solving system (1)-(3).

2. 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) were used. For the sake of simplicity, we use only square and cubic domains, specifically

and

.

2.1. Finite Difference Method

In the two-dimensional case, the domain is covered by a regular orthogonal mesh of nodeswith the step. Let us denote

(and 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

for;

This system is supplemented with the initial conditions

where.

In the case of the cubic domain, the spatial discretization is derived analogically. These two schemes will be called MOL.

2.2. Galerkin Finite Element Method

In this method, both the linear (triangular) and the bilinear (square) Lagrange elements (see [9]) are used. The basis of the finite dimensional subspace

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 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

(9)

where N_h is the number of the nodes of the spatial mesh, Equations (5)-(7) are transformed as follows:

(10)

(11)

(12)

for, where stands for the interpolant.

This system is supplemented with the initial conditions

where.

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-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.

3. 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

(13)

where the additional terms and assure that the system possesses the solution

where are positive constants.

Therefore, it turns out that

and (from the boundary conditions).

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

(14)

where, represents the numerical solution corresponding to the parameter 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:

(15)

where stands for a final time, and is a time step between the compared time levels of the solution. Two norms were used:

with

where denote the nodal values of the projections of and onto a regular orthogonal mesh with a spatial step.

The parameters were chosen as follows: L = 1, ,

, , , and. Further, (for the schemes FEMRHS, FEMHS and FEMSS) or (for the scheme MOL),; and the mesh for the evaluation of the error had nodes.

The results are summarized in Tables 1-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,

FEMSS appears to be the most accurate of all.

4. 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

and (similarly as in the articles [2] or [7]); and the initial conditions were of the form

The function was chosen constant. Therefore, the pattern types vary depending on the values of the parameters and.

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,). The value of the - norm 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.

5. Three-Dimensional Pattern Formation

This section provides examples of the spatial patterns that arise if we consider the cubic domain with or. Similar types of the initial conditions and values of the parameters as in the two-dimensional case were used, i.e., , and;

(denotes the third spatial coordinate). The values of and 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.

6. 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.

7. 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.

REFERENCES