An Efficient Explicit Finite-Difference Scheme for Simulating Coupled Biomass Growth on Nutritive Substrates

A novel explicit finite-difference (FD)method is presented to simulate the positive and bounded development process of amicrobial colony subjected to a substrate of nutrients, which is governed by a nonlinear parabolic partial differential equations (PDE) system. Our explicit FD scheme is uniquely designed in such a way that it transfers the nonlinear terms in the original PDE into discrete sets of linear ones in the algebraic equation system that can be solved very efficiently, while ensuring the stability and the boundedness of the solution. This is achieved through (1) a proper design of intertwined FD approximations for the diffusion function term in both time and spatial variations and (2) the control of the time-step through establishing theoretical stability criteria. A detailed theoretical stability analysis is conducted to reveal that our FD method is indeed stable. Our examples verified the fact that the numerical solution can be ensured nonnegative and bounded to simulate the actual physics. Numerical examples have also been presented to demonstrate the efficiency of the proposed scheme. The present scheme is applicable for solving similar systems of PDEs in the investigation of the dynamics of biological films.


Introduction
1.1.Background.The formation and dynamic evolution of biological films systems, such as bacterial biofilms, are generally very complicated.The investigation on the dynamics of biological films is essential for the development of biomedical industry, and it is becoming a rich and challenging topic of research in both the natural and the mathematical sciences.In order to better understand the mechanism, useful mathematical models have already been developed and are widely used to predict the dynamics of biological films systems [1][2][3][4][5][6][7][8].The features of the developing bacteria in the biofilm depend heavily on the environmental conditions, such as the relevant nutrition, temperature, moisture, and the oxygen level.If the substrate contains the nutrients that are beneficial to the colony, the development of biofilm will be promoted.On the other hand, the growth of the colony will consume the nutrients, which in turn decelerates the growth.Therefore, the biofilm development and the nutrients evolution form an elegant dynamic balance, known as the dynamics of biological films systems.To capture and predict the developing features of the biofilms systems, a lot of experimental studies have been carried out and many computer-based simulation techniques have been proposed.Mathematical models have also been developed to describe the pattern formation of photosynthetic bacteria on solid surfaces of flat-panel photobioreactors [9] to predict the photosynthetic bacteria (PSB) Rhodopseudomonas palustris CQK 01 biofilm formation [10], to study the survival and disinfection of Mycobacterium mucogenicum in potable water [7], and to propose the regulation of extracellular polymeric substances through quorum sensing in colonies of bacteria [8].Many of available mathematical models are expressed in terms of partial differential equations (PDEs) with a nonlinear density-dependent diffusion reaction describing the biofilm growth [9][10][11][12][13][14][15][16][17][18].A system of parabolic partial differential equations has been established to describe the complex interaction of a biomass system and a surrounding substrate of nutrients [12].Eberl et al. [1] proved the existence and uniqueness of nonnegative and bounded solutions of the problem, under suitable analytical constraints.
Inspired by Ebert's work, the present work aims to develop a novel explicit finite-difference (FD) method to effectively solve the nonlinear system of parabolic PDEs.A numerical scheme has also been proposed to solve the nonlinear PDE [2,5,6,12,18], revealing the sophisticated dynamics of biological films systems patterns.The present scheme is capable of producing an accurate solution approximating the positive and bounded growth of biological films.This is achieved through (1) a proper design of intertwined FD approximations for the diffusion function term in both time and spatial variations and (2) the control of the time-step through establishing theoretical stability criteria.A detailed theoretical stability analysis is conducted to reveal that our FD method is indeed stable.Our examples verified the fact that the numerical solution can be ensured nonnegative and bounded to simulate the actual physics.Furthermore, by choosing proper time-space grid, macro error can reach a certain precision.Numerical results of examples have also been presented to demonstrate the efficiency of the proposed scheme.The present scheme is applicable for solving similar systems of PDEs in the investigation of the dynamics of biological films.
Equation ( 1) is subjected to initial-boundary conditions given as follows: The long term behavior of the mass evolution of the initialboundary value problem associated with (1) was studied by others such as in [6].These useful studies have revealed the fact that for initial data  0 ,  0 satisfy the following conditions: Then, there exists a unique solution to (1)-(3) in a class of functions of the following: (1) ,  ∈  ∞ ( + × Ω) ∩ ([0, ∞),  2 (Ω)); (2) , () ∈  ∞ ( 1 (Ω),  + ) ∩ ( 2 (), [0, ∞)); (3) 0 ≤ (, ), (, ) ≤ 1, for every (, ) ∈ Ω ×  + and where which come from Eberl et al. [1] who used several mathematical assumptions observed experimental results.Theory evinces the fact that the numerical solution is nonnegative and bounded in the practice.In this paper, an efficient explicit FD method is introduced to approximate positive and bounded solutions of the dynamics of biological films systems with a nonlinear density-dependent diffusion reaction.We will show, through both theoretical analysis and numerical example, that the solution process is effective, and the positivity is always.
The rest of the presentation is organized as follows.We introduce two explicit positivity-preserving finite-difference schemes in Section 2. The properties of these two schemes will be examined in Section 3. In Section 4 four example problems are simulated using the present method, and the computational efficiency is demonstrated.Finally, we conclude with some remarks in Section 5.

Explicit Positivity-Preserving FD Scheme: A Novel Technique
We introduce now an explicit positivity-preserving FD scheme for the problem defined in (1)- (3).Consider a regular space-time domain: where governing equations set (1) can be expressed in the following form: We first note that the first-order difference approximation can be written in a number of ways using a regular grid defined in the space-time domain with intervals of Δ, Δ, and Δ.For example, the time derivative can be given in a standard FD form as The difficulty is how to make use of these approximations to the nonlinear PDEs (6), so that it can be linearized effectively and produces stable and convergent numerical solutions.Our idea is to explore these various FD schemes that work best for the unique structure of our PDE.We use first the following approximation for the diffusion function: () ≈  ( ( ± Δ, , )) −  ( (, , ))  ( ± Δ, , ) −  (, , ) .
To get rid of the denominator in the above equation and simplify our discrete equation system, we now choose to use the following approximation for the square of the first derivative: ( ± Δ, , ) −  (, , ) Δ  ( ± Δ, , ) −  (, ,  + Δ) Δ . ( Combining these approximations given in ( 9) and ( 10),   ()(  ) 2 can now be given as It is clear that the complexity of the nonlinearity is significantly reduced to a manageable form.By the same token, we will have To further simplify the notations, ( 11) and ( 12) can be denoted, respectively, using  ±  and  ±  , and we can denote ((, , )) by   .Equations set (6) can now be rewritten in a concise form of Let , , , and  be real numbers such that  <  and  < .An open rectangle problem domain can be defined as Ω = (, ) × (, ) in  2 .A set of regular grids in the spatialtime domains [, ], [, ], and [0, ] can be easily generated as where , , , , , and  are all positive integers.
These constraints may be homogeneous Neumann or Dirichlet conditions and can be expressed through the following unified discrete forms: If the constant  is equal to 1, then the homogeneous Neumann boundary condition was imposed on Ω; if  is set to 0, then the homogeneous Dirichlet boundary condition is imposed on Ω.
The discrete form of initial conditions will be expressed through the following forms: In order to estimate the accuracy of the presented scheme, we use the following error estimation to measure the error of numerical solution.
For the homogeneous Neumann boundary condition, we integrate equation system (6) on the domain (, , ) ∈ [0, ] × [0, ] × [0, ].We have the following forms: Letting () = ∫  0 (∬ Ω (/( 4 + ))Ω) in equation system (18), eliminating the term (), we have To further simplify the notations, let Thus, (19) can be rewritten as the following form: The macro error can be defined as where The term   (, , ) is the numerical solution of substrate concentration containing the nutrients and term   (, , ) is the numerical solution of density of biomass.

Properties of the Scheme: A Theoretical Analysis
In this section, we give sufficient conditions which ensure that the numerical solutions of our schemes are nonnegative and bounded.For the simplicity in notation, let Thus, ( 15) can be rewritten as in the following form: Lemma 1 (nonnegativity and boundedness of Scheme ( 15)). ( In the following, we prove that  +1 , is nonnegative.Proof. +1 , < 1 is clear. +1 , ≥ 0 is equivalent to the right-hand side of Scheme (15), and hence it is greater than or equal to zero; that is, The above can be rewritten in the following form: Proof. +1 , ≥ 0 is clear. +1 , < 1 is equivalent to the right-hand side of Scheme (15), and hence it is less than 1; that is, The above can be rewritten in the following form: Therefore, (i) if ( 2 − 2 )  , >  2 holds, the above is equivalent to Thus, if   , < 1 and the conditions of (i) hold, the numerical solution according to Scheme (15) is always less than 1.
Remark 3. In Scheme (15), if the decay rate of biomass  2 is greater than the maximum specific growth rate  3 , then the biomass  is in attenuation governed by a nonlinear advection-diffusion equation.Thus, our Scheme ( 15) is unconditionally stable.We conclude that our method can always ensure stability and boundedness independent of the choice of the space-time grid.On the contrary, if the decay rate of biomass  2 is less than the recombination coefficient  3 /( 4 + ), then the biomass  is growth that is governed by a nonlinear advection-diffusion equation.

Numerical Examinations
We now provide some numerical examples to examine the performance of our method.Consider a nondimensional spatial domain Ω = [0, 1] × [0, 1].Our numerical tests will be conducted for biomass growth on nutritive substrates with several of initial density profiles in the domain.
The detailed dynamic evolution of the biomass will be computed as a history of time.Four typical examples are studied: (1) coupled biological density attenuation system; (2) coupled biological density grown system; (3) complex coupling system, in which the density of biological evolution process is increasing at the first stage and then decreases; (4) coupled biological system, in which discontinuous initial conditions are imposed, with homogeneous Dirichlet boundary conditions on the all boundaries.
In Figure 1, it can be observed that the biofilm density continuously decays and the colony spreads with the increase of time.This is because the biomass decay rate  2 is larger than the maximum specific growth rate  3 .We have proved that the case is an attenuation model.By comparing with the distribution graph of biofilm density and substrate concentration of nutrition at time levels  = 3, 7, and 10, a complementary phenomenon can be observed; that is, the density biofilm is relatively large in the domain where the nutrient consumption is relatively large; the density of biofilm is relatively small in the domain where the nutrient consumption is relatively small.Furthermore, the figure on the left and the right figure are highly complementary to each other.It confirms that the development process of biofilm density and nutrient density is highly coupled.
Table 1 lists the macro error at three different times for different .From Table 1 and Figure 2, the following can be found.
(1) When  is less than or equal to 1, macro error can reach to 10 −3 for Example 1. (2)  is smaller and macro error is smaller.This means that our model is convergent for variable .
Example 2. Let Ω = (0, 1)×(0, 1) and  = (, ).We consider an initial density profile with the form of where  1 = 0.47,  2 =  4 = 0.5, and  The results are plotted in Figure 3, using Δ = 0.00004.It is seen that the biofilm density grows continuously and spreads with the increase of the time.This is because the biomass decay rate  2 is much smaller than the maximum specific growth rate  3 .The biofilm density will attenuate at later time, which will be further examined in Example 3.
In comparing the density distribution of the biomass with substrate concentration of nutrition at time levels  = 3, 7, and 10, the same complementary phenomenon can also be observed; that is, the biomass density is relatively large at in the domain where the nutrient consumption is relatively large and relatively small in the domain where the nutrient consumption is relatively small.Comparing the figure on the left with that on the right, it is seen that the development process of biofilm density and nutrient density is highly coupled.Table 2 lists the macro error at three different times for different .From Table 2 and Figure 4, the following can be found.
(1) Macro error can reach to 10 −3 at  = 3, 7 and reach to 10 −2 at  = 10 for Example 2. (2)  is smaller and macro error is smaller.This means that our model is convergent for variable .
From Figure 5, it can be seen that the numerical solution which is obtained using the traditional Euler difference method is unstable, and the numerical solution becomes unbounded at time levels  = 7, 10, but presented Scheme (15) can ensure that the numerical solution is nonnegative and bounded.Therefore, our Scheme ( 15) is stable for Example 2.
From Figure 6, we observe that the biofilm density grows fast and spreads over time from its initial profile.It then starts to decay on some domain.This is because the biomass decay rate  2 is smaller than the combined effective coefficient  2 =  3 (  , /( 4 +   , )).Because the nutrition function (, , ) is attenuating, which brings about the decrease of the recombination coefficient  2 =  3 (  , /( 4 +   , )) over time, the biomass decay rate  2 will be larger than the combined effective coefficient  2 =  3 (  , /( 4 +   , )).At later time, the biofilm density will shift from the growth process to the attenuation process.By comparing the density distribution graph of biomass with substrate concentration of nutrition at time levels  = 4, 8, and 12, the development of biofilm density and nutrient density has a high degree of correlation, which is consistent with Examples 1 and 2. The figures in the two columns are mutually complementary, which further shows that our model is reasonable, reliable, and effective.
Table 3 lists the macro error at three different times for different .From Table 3 and Figure 7, the following can be found.
(1) Macro error can reach to 10 −3 at  = 4, 12 and reach to 10 −2 at  = 8 for Example 3. (2)  is smaller and macro error is smaller.This means that our model is convergent for variable .
From Figure 8, it can be observed that the numerical solution which is obtained using the traditional Euler difference method is unstable, and the numerical solution becomes unbounded at time levels  = 8, 12, but presented Scheme (15) can ensure that the numerical solutions are nonnegative and bounded.Therefore, our Scheme ( 15) is also stable for Example 3.
Oscillations are not observed in present Scheme (15) but clearly observed in standard EE method (seeing it in Figure 10).
The above numerical simulations give four different types of development process of biofilm density and nutrition density.We capture the evolution and cloning features of biological films dynamics systems, which are consistent with the actual development process of biological films dynamics systems.

Conclusions
In this work, we present a novel explicit finite-difference method to approximate the positive and bounded solution of coupled substrate-biomass system, which is governed by two parabolic partial differential equations with a nonlinear density-dependent diffusion reaction and two unknown functions.Theoretically, the existence and uniqueness of nonnegative and bounded solutions of the model have been proved in the specialized literature, while acquiring the exact analytical solution of the biological model is very difficult.Therefore, we have to resort to approximate solutions using numerical simulation methods.This paper presents a welldesigned explicit finite-difference scheme for this very purpose.Our approach in designing such a scheme gives a new idea on how to transfer nonlinear terms in the PDE into linear ones in smart and effective ways, while ensuring the stability of the solution and the bounded properties.Our theoretical and numerical studies have concluded the following.
(1) In the practical applications, we should choose proper time-step for desired accuracy in numerical solutions.The determination of the "proper" spatial and time grids depends on the actual problem and needs to be done via essential trial and error.
(2) All these examples have showed that our novel methods can always ensure the solution is automatically bounded within zero and one, and the development process of biofilm density and nutrient density is highly coupled.
(3) Our presented Scheme can eliminate the phenomenon of oscillation for the initial value problems of discontinuous.

Figure 2 :
Figure 2: Macro errors versus different  at three different times.

Figure 3
Figure 3 plots the distribution graph of biological and substrate concentration of nutrition at four different times.It is seen that the biofilm density grows continuously and spreads with the increase of the time.This is because the biomass decay rate  2 is much smaller than the maximum specific growth rate  3 .The biofilm density will attenuate at later time, which will be further examined in Example 3.In comparing the density distribution of the biomass with substrate concentration of nutrition at time levels  = 3, 7, and 10, the same complementary phenomenon can also be observed; that is, the biomass density is relatively large at in the domain where the nutrient consumption is relatively large and relatively small in the domain where the nutrient consumption is relatively small.Comparing the figure on the left with that on the right, it is seen that the development process of biofilm density and nutrient density is highly coupled.

Figure 4 :
Figure 4: Macro errors versus different  at three different times.

Figure 7 :
Figure 7: Macro errors versus different  at three different times.
1.2.Briefing on the Mathematical Model.The positive and bounded evolution of biomass and substrate of nutrients on biological films are governed by the following coupled system of PDEs: All these coefficients are obtained by experiment.The term () represents the diffusion coefficient of biomass, that is, in general, a function of biomass density, and hence it is a source of the nonlinearity.Experimental observation has found that the diffusion coefficient degenerates for small biomass densities and is singular at the upper bound of the biomass density.The form of the diffusion coefficient is given as a function of the biomass itself: 1 ∇ 2  (, ) −  1  (, )  (, )  4 +  (, ) ,   (, ) =  2 ∇ ⋅ ( ( (, )) ∇ (, )) −  2  (, ) +  3  (, )  (, )  4 +  (, ) .(1) Here, the spatial operators ∇ and ∇ 2 denote, respectively, the gradient and the Laplace operators, where Ω in   ( = 1, 2, 3) is an open domain, where the biomass grows, Ω is the boundary of Ω, and Ω = Ω ∪ Ω is the closed domain, where the biomass growth problem is set.(, ) is the unknown density of biomass that is normalized with respect to the maximum biomass density over the entire space-time domain Ω×[0,  max ], (, ) is the unknown 2 : biomass diffusion coefficient;  1 : maximum specific consumption rate;  2 : biomass decay rate;  3 : maximum specific growth rate;  4 : Monod half saturation constant.