Numerical approximation of one model of bacterial self-organization

Abstract. This paper presents finite difference approximations of one dimensional in space mathematical model of a bacterial self-organization. The dynamics of such nonlinear systems can lead to formation of complicated solution patterns. In this paper we show that this chemotaxisdriven instability can be connected to the ill-posed problem defined by the backward in time diffusion process. The method of lines is used to construct robust numerical approximations. At the first step we approximate spatial derivatives in the PDE by applying approximations targeted for special physical processes described by differential equations. The obtained system of ODE is split into a system describing separately fast and slow physical processes and different implicit and explicit numerical solvers are constructed for each subproblem. Results of numerical experiments are presented and convergence of finite difference schemes is investigated.


Introduction
Many mathematical problems of biological systems are described by non-stationary and non-linear diffusion-advection-reaction equations.The dynamics of their solutions can be very complicated, the interaction of different physical processes can lead to development of spatial and temporal patterns and instabilities [1].The development of complicated spatial-time patterns was observed also in real bioluminescence images (see [2,3]).
In addition, solutions of biological and chemical systems satisfy special properties, such as positivity, boundedness and conservativity [4].Therefore the development of robust and efficient numerical algorithms for solution of such problems still remains a very important challenge of computational mathematics.
A substantial number of research exists on systems related to the chemotaxis-growth systems.Here we comment on some specific results, closely connected to our analysis (see the paper of Painter and Hillen [5] for a recent excellent review on this topic).
It is well-known that the theory of exponential attractors explains important properties of dynamical systems in infinite-dimensional spaces.The exponential attractor is c Vilnius University, 2012 a compact set with finite fractal dimension which contains a global attractor interiorly and attracts every trajectory in an exponential rate.Exponential attractors are also known to have very strong stability in approximation.This gives a possibility to show a global reliability of numerical computations.Recently very important results on the fractal dimension estimate of the global attractors for abstract quasilinear parabolic evolution equations are obtained by Aida et al. [6].Application of these methods for chemotaxis-growth systems is done by Nakaguchi and Efendiev [7].Using non-negativity of solutions they managed significantly to improve dimension estimates with respect to the chemotactic parameter.
Many models of chemotaxis have been formulated [8].Here we note the paper of Painter and Hillen [5], where the dynamic properties of one-dimensional Keller and Segel model are deeply investigated.Their analysis gives a template for exploration of different models of chemotaxis: the linear stability of the homogeneous steady state is investigated, then it is shown that most of multiple-peak patterns are unstable and a coarsening process leads to a unique global aggregation.It is revealed that the long-time dynamics of the solutions fall into on of the four classes: homogeneous steady state, stationary spatial patterns, spatio-temporal periodic and irregular solutions.We note also the analysis of numerical stability for the various solution classes, when the influence of small perturbation applied at some time moment is computed by comparing perturbed and unperturbed solutions.The numerical Lyaponov exponent is used to get a quantative measure for this analysis.
It is well-known that solutions of the one-dimensional case may diverge in a finite time only if the chemoattractant does not diffuse, but in higher dimensional spaces divergences are common.The Keller and Segel model allows pattern formation and spatiotemporal chaos in one spatial dimension, but these features have not been investigated in full detail even in the (2+1)-dimensional case [9].Thus the development of efficient numerical algorithms (including parallel algorithms) for simulation of multidimensional chemotaxis models is an important task.Efficient and accurate finite element solvers for two-dimensional chemotaxis problems are proposed and investigated by Strechl et al. [10].They have analyzed the monolithic and decoupled variants of different algorithms and applied various techniques for linearization of the obtained discrete systems.
In this paper we consider the development of numerical algorithms to solve one system of PDE proposed in [11] for simulation of the bacterial self-organization in circular container.We suggest to explain a chemotaxis-driven instability by connecting it to the ill-posed backward in time diffusion process.The main goal of this paper is to construct robust and efficient numerical approximations.These algorithms are targeted for development of parallel algorithms in the case of 2D and 3D problems.The second goal is to to investigate the convergence of proposed finite-difference schemes in the case of irregular dynamical behaviour of solutions.
The rest of the paper is organized as follows.In Section 2 the mathematical model for simulation of the bacterial self-organization in circular container along the contact line is presented.Linear stability analysis of the model is done and it is shown that chemotaxisdriven instability can be connected to the ill-posed problem defined by the backward in time diffusion process.In Section 3 we give details of the numerical techniques which are used for the construction of finite-difference schemes.Apriori estimates of the boundedness and positivity of the discrete solution are proved.Results of numerical experiments are presented in Section 4. The convergence of the discrete solution is investigated in the case of regular and irregular dynamical behaviour of the nonlinear system.For irregular case, it is proposed to consider shorter time intervals for the classical convergence analysis, then the accuracy of the finite difference scheme can be investigated by using the Runge rule.Finally, in Section 5 some conclusions are formulated.

Problem formulation 2.1 Mathematical model
We consider a mathematical model for simulation of the bacterial self-organization in circular container along the contact line, which is proposed in [11].An excellent review on PDE models for chemotaxis can be found in [8] and references contained therein.
By using some simplifying assumptions, the mathematical model is defined as a system of one dimensional nonlinear equations: where u is the dimensionless cell density, v is the dimensionless chemoattractant concentration, α ≥ 0 defines the receptor sensitivity, β > 0 stands for saturating of the signal production, γ > 0 defines the ratio of the spatial and temporal scales, and p ≥ 1.Thus a simple model with the linear production of chemoattractant is replaced with a saturating term which depends on u in nonlinear way, and a chemotactic sensitivity coefficient is also modified to take into account a quadratic saturation with respect to v.
The initial conditions are defined as The periodicity conditions are formulated as boundary conditions: (3)

Analysis of the mathematical model
In [11], the non-deterministic initial conditions are investigated: where ε(x) is a 20% random uniform spatial perturbation.
Nonlinear Anal.Model.Control, 2012, Vol.17, No. 3, 253-270 A random perturbation of initial data is not a necessary condition for self-organization of the solution.The sensitivity of the solution with respect to initial data should be investigated, i.e. well-posedness of the mathematical model should be analysed.Many technological and physical processes can lead to development of spatial and temporal instabilities in solutions, we mention only a fingering instability in thin evaporating liquid films [12,13], in buoyancy-driven fluid filled cracks [14] or in porous media flows [15], and pattern formation in reaction-diffusion systems (see a fundamental book of Murray [16] and references given in it).The chemotaxis-driven instability is investigated in many papers and different types of dynamic behaviour of nonlinear systems can be observed, see, e.g., the classical Keller-Segel model [16,17].
The dynamics of nonlinear systems can be investigated by various mathematical techniques, an extensive review with many examples and applications can be found, e.g., in [1,18] and references contained in these books.

Linear stability analysis
Here we apply a standard linear stability analysis around the homogeneous steady state solution [5,16].The steady state of model ( 1) is given by ū = 1, v = 1/(1 + β).Linearization around this solution gives the following system: for small perturbations U (x, t), V (x, t).The stability of the steady state depends on the temporal eigenvalues of the stability matrix where k ≥ 0 denotes the wavenumber.In the case of zero-flux boundary conditions and the interval [0, L] we have k = nπ/L, n ≥ 0. If the matrix A k has at least one eigenvalue with a positive real part, the homogeneous steady state is unstable [5,16].The eigenvalues λ of A k are obtained from the characteristic equation Rescaling k 2 = γl 2 , we can use results of stability analysis in [5].Then we obtain the following necessary condition for stability matrix A k to have an eigenvalue with a positive real part: We see that this condition does not depend on parameter γ.

Backward in time diffusion
In this paper we apply a new approach, when the chemotaxis-driven instability is correlated to the ill-posed problem defined by the backward in time diffusion process.In general, the well-posedness of the model is connected to the important property of chemotaxis process, that the velocity of advection of u depends on the gradient of the chemoattractor.Thus, if the chemo-attractor depends monotonically on the cell density, then such a dependence leads to anti-diffusion flow of the cell density.For simplicity of analysis, let us assume that parameter γ is sufficiently large, but γr ∼ O(1), then due to the fast relaxation we get that Thus in the limit case, equation ( 1) can be written as a nonlinear diffusion-reaction equation For some values of parameters equation ( 4) describes a backward in time parabolic problem and therefore it is ill-posed [19,20].Solution of ( 4) does not depend continuously on initial data and small perturbations can lead to unbounded changes of the solution.Such a scenario potentially can describe the dynamics of the model (1).We note that values of parameters, which were used in [11] for computational simulations, belong to this critical set.
Next we present a simple analysis of the backward parabolic equation.As an example we consider a linear backward parabolic equation with the Dirichlet boundary conditions Let define the initial condition u 0 (x) in ( 5) as a solution of the forward parabolic equation: at t = T , i.e. u 0 (x) = v(x, T ).The well-known properties of the forward parabolic equation are straightforwardly obtained by applying the standard Fourier analysis.Let assume that the initial condition can be written as v(x, 0) = N n=1 c n sin(πnx).Then a solution of parabolic problem ( 6) is given by Nonlinear Anal.Model.Control, 2012, Vol.17, No. 3,[253][254][255][256][257][258][259][260][261][262][263][264][265][266][267][268][269][270] We see that high modes of the Fourier sum are damped faster and v(x, t) becomes smoother for t > 0. Fig. 1 shows the initial distribution v(x, 0) and the solution v(x, T ) of classical parabolic problem (6) at time T = 0.1.The diffusion coefficient D = 0.1 is used in all computations.The Fourier stability analysis of the backward parabolic problem can be done in a similar way: Now we see that high modes of the Fourier sum are growing faster and noise perturbations (e.g., corresponding to a white noise) are amplified for t > 0.
Fig. 2(a) shows a solution u(x, 0) of the backward parabolic problem at time moment t = 0 (the initial condition u(x, T ) is defined by the solution v(x, T ) of problem ( 6)).All useful information is lost due to amplification of the noise modes.
Regularization of the backward parabolic equation can be done in various ways.Here we mention two approaches.In the first method the regularization is done by a nonlocal boundary value problem [21]: with a > 1 being given data and α > 0, the regularization parameter.As an approximation of u(x, t) we take u α ((a − 1)T + t).A priori and a posteriori rules for a selection of regularization parameter α are proposed in [21], which yield order optimal regularization methods.The standard Fourier analysis gives an explicit form of the solution from which the well-posedness of the regularized problem (7) follows. www.mii.lt/NA In the second approach the ill-posed backward parabolic problem is regularized by a perturbed well-posed backward parabolic equation [20]: with α > 0 being the regularization parameter.The Fourier analysis explains the regularization effect: We see that high modes of the Fourier sum are damped fast and therefore noise perturbations are not amplified for t < T .Fig. 2(b) shows a solution of the regularized problem (8).
We note, that due to non-linearity of the diffusion coefficient, the ill-posedness of the problem is self-limiting, since 3 The finite difference scheme In this section we present numerical techniques which are used to approximate solutions of system (1)- (3).A comprehensive treatment of theoretical and implementation issues of discretization methods for advection-diffusion-reaction problems are given in the monograph by Hundsdorfer and Verwer [4].Very interesting applications of these results for biomedical problems are described by Gerisch and Chaplain [22], Gerisch and Verwer [23] (see also references given in these publications).Efficient and accurate finite element solvers for chemotaxis problems are proposed and investigated by Strechl et al. [10].They have compared the monolithic approach and the decoupled variant of algorithms.In particular they have analyzed stability, accuracy and efficiency of the different methods.We note that such results are important when a parallelization strategy should be selected.

The method of lines. Discretization in space
We use the method of lines (MOL) approach (see [4,22]).At the first step we approximate the spatial derivatives in the PDE by applying robust and accurate approximations targeted for special physical processes described by differential equations.Domain [0, 1] is covered by a discrete uniform grid with the grid points x j .On the semidiscrete domain ω h (k) × [0, T ] we define functions U j (t) = U (x j , t), V j (t) = V (x j , t), j = 0, . . ., N −1, here U j , V j approximate solutions u(x j , t), v(x j , t) on the discrete grid ω h at time moment t.We also define the forward and backward space finite differences with respect to x: Using the finite volume approach, we approximate the diffusion and reaction terms in the PDE system (1) by the following finite difference equations: The stencil of the discrete scheme requires to use functions defined outside of the grid ω h .We apply periodicity boundary conditions (3) to define discrete functions for j < 0 or j ≥ N : The advection term in equation ( 1) depends on the variable velocity therefore the maximum principle is not valid for the respective transport equation.But problem (1) still has non-negative solutions, and this property can be preserved on the discrete level by applying proper upwinding approximations.The discrete spatial approximation of the velocity is computed by In the following we consider the upwind-based discrete fluxes [4,22]: with the Koren limiter function The limiter ψ depends on the smoothness monitor function For ψ = 0 we get the standard first-order upwind flux F TUW U, a, j + 1 2 = max a j+1/2 , 0 U j+1 + min a j+1/2 , 0 U j .

Let us denote the discrete advection operator as
Then we obtain a nonlinear ODE system for the evaluation of the approximate semidiscrete solutions

Operator splitting methods
In order to develop efficient solvers in time for the obtained large ODE systems we apply the splitting techniques.They take into account the different nature of the discrete operators defining the advection A T (U, V, j) and the diffusion-reaction A DR1 (U, j), A DR2 (V, U, j) terms.The system resolving the semi-discrete advection process can be solved very efficiently by using explicit solvers, while the diffusion-reaction semi-discrete system is stiff and it requires an implicit treatment.Also we are interested in preserving at the discrete level the positivity and/or boundedness of the solution, if such properties hold for the differential ODE system.The splitting method gives us a possibility to construct robust and flexible parallel algorithms, when domain decomposition should be incorporated into numerical approximations.First we consider the symmetrical splitting method (also known as the Strang splitting [24]).Given approximations U n j , V n j at time t n , we compute solutions at t n+1 = t n + τ by the following scheme: Here we split the given ODE system into two blocks with respect to U j and V j functions.
Lemma 1. Solutions of the splitting ODE problem (9)-( 12) are non-negative if U n j ≥ 0 and V n j ≥ 0 for all x j ∈ ω h .Proof.The proof for the advection subsystems follows from the construction of the discrete fluxes by using the upwinding technique.The proof for the diffusion-reaction subsystems follows from the lemma in [23] that the solution of an initial value problem for systems of ODEs is non-negative if and only if for all t and any vector V ∈ R m and all It is easy to see that for the diffusion-reaction subsystems the diffusion parts of the matrices are diagonally dominant and all off-diagonal entries are non-positive.For the reaction parts the required estimates are also trivially satisfied.
Lemma 2. If 0 ≤ U n j ≤ C for all x j ∈ ω h , then a solution of the splitting ODE problem (10) is also bounded U j ≤ max(C, 1).
Proof.Here we use the fact that U = 1 is a stable attractor of the reaction function.Let U n i = max j U n j .Then it follows from the definition of A DR1 (U, j) that in the worst case The lemma is proved. www.mii.lt/NA

Numerical integration of ODEs
There are many numerical integration methods for solution of non-stiff and stiff ODEs.
For detailed discussions of these schemes we refer the reader to [4,22,[25][26][27].Let ω τ be a uniform time grid here τ is the time step.For simplicity this step size is taken constant.
In the following, we consider numerical approximations U n j , V n j to the exact solution values u(x j , t n ), v(x j , t n ) at the grid points (x j , t n ) ∈ ω h × ω τ .
Remark 1.In [11], the explicit forward Euler scheme is used to solve problem (1).Since no details are given in [11] on approximations of spatial derivatives, we use discrete operators introduced in previous sections and write the explicit forward Euler scheme as: We note that this scheme can be written as a splitting algorithm: Thus the explicit Euler scheme can be considered as a special case of splitting algorithms.Despite easy implementation and good parallel properties of explicit algorithms, the main drawback of the explicit Euler method is that due to the conditional stability we must restrict the integration step to τ ≤ Ch 2 for stiff discrete diffusion-reaction subproblems.
Here we propose to use a linearized implicit backward Euler method for the approximation of the diffusion-reaction subproblems and the explicit forward Euler method for solution of the advection subproblem.We have restricted to the first order methods due to their robust stability.Note that our main goal is to investigate the influence of a possible ill-posedness of the PDE system to the asymptotical behaviour of the solution.
We discretize the semidiscrete problem ( 9)-( 12) with the fully discrete scheme We apply two splittings of the advection term, because then we use only half of the splitting step size for the explicit method.This doubles the stability and positivity domains of the explicit method (see [23]).
Next we prove that statements of Lemma 1 and 2 hold also for solutions of the fully discrete finite difference scheme ( 13)-( 16).
Lemma 3.For a sufficiently small time step τ ≤ τ 0 solutions of the finite difference scheme (13)-( 16) are non-negative if U n j ≥ 0 and V n j ≥ 0 for all x j ∈ ω h .Proof.The proof for the advection problems ( 13) and ( 15) follows from the construction of the discrete fluxes by using the upwinding technique and selection of a sufficiently small time step τ ≤ τ 0 .
The proof for the diffusion-reaction problems ( 14) and ( 16) follows from the maximum principle [28].For example, consider problem (14).We assume, that We write the discrete equation ( 14) for U n+2/3 i in an explicit form ≤ C for all x j ∈ ω h , then a solution of the finite difference scheme ( 14) is also bounded Proof.The proof is based on the maximum principle and a special form of the discrete reaction term.First, we consider the case www.mii.lt/NANext we consider the case C > 1.Then we get that The lemma is proved.
We note that a convergence analysis of linear splitting schemes is quite well-developed for many classes of physical processes, see [4,28].For nonlinear problems such results are proved only for some partial cases.We will deal with the convergence questions in a separate paper.

Numerical experiments
In this section, we present results of numerical experiments in order to verify our theoretical investigations.

Convergence analysis
The first goal is to investigate the convergence of the finite difference discretizations of the PDE model (described in the previous section) in the case when dynamics of solutions leads to formation of complicated spatial-temporal patterns.The development of the bacterial population was simulated for the following values of the model parameters (see [11]): First, we consider the MOL type discrete scheme ( 13)-( 16), when the space grid step has been fixed to h = 0.005, and computations has been performed with different time steps τ = 1•10 −6 , 5•10 −7 , 2.5•10 −7 .A simple initial condition u 0 (x) = 1+0.2sin(4πx), x ∈ [0, 1] is used in computations.Fig. 3(a) shows the snapshots of function U j (t n ) at time T f = 0.6.
As can be seen, the discrete solutions do not converge when the discretization temporal step τ is decreased.
A similar behaviour of numerical approximations is also observed for the fully discrete scheme ( 13)- (16).Fig. 3(b) shows the results, when the spatial and temporal grid steps are connected by the relation τ = 5 • 10 −5 h and h = 0.01, 0.005, 0.0025.
In order to prove that such convergence behaviour of discrete solutions depends on dynamical instability of the mathematical model ( 1) and more exactly it is a consequence of negative diffusion driven instability due to chemotaxis process, we consider the same problem (1) but for a short time interval t ∈ (0, T f ].During this time the diffusion coefficient D still dominates the chemotaxis perturbations and the mathematical model defines a well-posed problem.Dynamics in time of the maximal advection velocity v max is presented in Fig. 4(a), thus we can take T f = 0.125.Next we have investigated a sensitivity of the numerical solution with respect the initial conditions computed at some time moment, when the solution have reached a phase of pattern formation.The main goal is to determine a length of time interval for which the discrete solution still remembers the influence of initial conditions.To address this, we have computed a solution of ( 13)-( 16) till t = 1.5 with H = 1/4000 and used it as initial conditions for the following computations.The error is defined as the difference of two solutions U h , U H computed with space steps h and H respectively: Table 1 presents the dynamics of errors Z h (T ) for different final times T and temporal and space mesh steps.Table 1.Sensitivity and convergence analysis of the discrete solution with respect to the initial condition.Errors of the solution of ( 13)-( 16) are presented for different temporal and space mesh steps at T = 1.52, 1.54, 1.58, 1.65.Discrete solution of ( 13)-( 16) at t = 1.5 with h = 1/4000 is used as the initial condition.The obtained results indicate that the proposed finite difference scheme is robust and accurate.Till time moment T = 1.58 the discrete solution converges when temporal and space mesh steps are reduced, the convergence order coincides with theoretical estimates.
We note that in [5] the numerical stability for various solution classes was investigated by considering the impact of a small perturbation applied at t = T p .The subsequent difference between perturbed and unperturbed solutions was tracked and analyzed.The numerical Lyaponov exponent was used to estimate the stability of the solution.

Formation of complex patterns
In order to study the dependence to initial conditions of pattern formation we have simulated the model ( 1)-( 3) with values of parameters (17).It is expected to get spatiotemporal irregular solutions and according to Aida et al. [6] the qualitative form of spatiotemporal patterns should not depend on specific initial conditions (see also [7]).Fig. 5 plots typical results for two types of initial conditions: u(x, 0) = 1 + 0.2 sin(2πx) and u(x, 0) = 1 + 0.2 sin(4πx).
The space and time dynamics of cell density is presented.
We see that a similar qualitative pattern formation is observed for different initial conditions.
In order to test the regularity of the solution, the corresponding phase plane trajectory (U (x, t), V (x, t) at x = L/2 is computed and plotted in Fig. 6 for initial condition u(x, 0) = 1 + 0.2 sin(2πx).It has a typical strange attractor appearance.
Finally, we have added the regularization term into the first equation of system (1).The simulations show that if the regularization parameter α is such, that the perturbed Nonlinear Anal.Model.Control, 2012, Vol.17, No. 3, 253-270 problem starts to be well-posed, the pattern formation disappears and the cell density converges to a homogeneous stable stationary state.(17).The time trajectory in the (U (x, t), V (x, t)) phase plane at x = L/2.It has a strange attractor appearance.

Conclusions
In this paper we have studied one-dimensional diffusion-chemotaxis-reaction mathematical model and have shown that chemotaxis-driven instability can be connected to the ill-posed problem defined by the backward in time diffusion process.By using the well-established techniques the mathematical model is approximated by the discrete computational model.It is proved that the discrete solution inherits main properties of the solution of the differential mathematical model.Results of numerical experiments show that for such problems, where chemotaxis-driven instability defines the dynamics of a solution, the classical convergence property of numerical algorithms is not applicable for long time intervals.We propose to make the classical convergence analysis in short time intervals, then the accuracy of the discrete schemes can be evaluated by the Runge rule.
Instead of pointwise and similar convergence metrics qualitative criteria or averaged statistical characteristics can be used, as in the discrete element method.
The one-dimensional model can be too simple in order to obtain a satisfactory agreement with pattern formation seen in experiments.Thus 2D and 3D generalizations of the given mathematical model should be investigated.Construction of robust and efficient parallel algorithms can be done by using domain decomposition and splitting in space methods [29].

Fig. 1 .
Fig.1.Plots of the initial condition v(x, 0) (black) and the solution v(x, T ) (red) of classical parabolic problem(6).Function v(x, T ) is used as initial condition for the backward parabolic problem(5).(Online version in colour.)

Fig. 2 .
Fig.2.Plots of solutions of the backward parabolic problem (5): (a) the initial state at t = 0.1 (red), the exact solution at the final time moment t = 0 (black) and the regularized solution of (8) at t = 0 (blue) when regularization parameter α is too small; (b) the same information for regularization parameter α = 0.9 • 10 −5 .All useful information is lost in the first case, while v(x, 0) is reconstructed well in the second case.(Online version in colour.)