#### Abstract

A novel explicit finite-difference (FD) method is presented to simulate the positive and bounded development process of a microbial 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.

#### 1. 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–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–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 a proper design of intertwined FD approximations for the diffusion function term in both time and spatial variations and 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.

##### 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: Here, the spatial operators and denote, respectively, the gradient and the Laplace operators, where in 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 is the unknown corresponding substrate concentration containing the nutrients, and , , , , , and are given constants with the following meanings: : substrate diffusion coefficient; : biomass diffusion coefficient; : maximum specific consumption rate; : biomass decay rate; : maximum specific growth rate; : Monod half saturation constant.

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: where and .

Equation (1) is subjected to initial-boundary conditions given as follows: The long term behavior of the mass evolution of the initial-boundary value problem associated with (1) was studied by others such as in [6]. These useful studies have revealed the fact that for initial data satisfy the following conditions:(a) and for every ;(b) and ;(c), for every and .Then, there exists a unique solution to (1)–(3) in a class of functions of the following:(1);(2);(3), 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.

#### 2. 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 At the current time level and new time level , the second-order spatial partial derivatives and can also be evaluated with the following forms of FD approximations: 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), 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 . A set of regular grids in the spatial-time domains , and can be easily generated as where , and are all positive integers.

When , let and . Let , , respectively, be the approximation to , and expresses the value of . The explicit positivity-preserving FD scheme can be designed as follows: for every , , and , where , .

It is now necessary to impose appropriate, discrete constraints on the edges of the rectangle domain . This can be expressed through the following forms.

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 . We have the following forms: Letting 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.

#### 3. 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)). *(1) If and hold, then the numerical solution is nonnegative and bounded. That is, , if and .**(2) If holds, then an arbitrary time-step can always insure , if and .*

*Proof. * is clear.

Hence,

In the following, we prove that is nonnegative.

is equivalent to ; that is, .

So, If and hold, then the numerical solution is nonnegative. If holds, then an arbitrary time-step can always insure , if and .

Lemma 2 (nonnegativity and boundedness of Scheme (15)). *(1) In this case, .*(i)*If and hold, then the numerical solution satisfies , if and .*(ii)*If holds, then an arbitrary time-step can always insure , if and .*

*Proof. * is clear.

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:
Therefore, (i) If and hold, then the numerical solution satisfies , if and .

(ii) If holds, then an arbitrary time-step can always insure , if and .

The limitation for time-steps set can insure that the numerical solution of biomass is nonnegative and such limitation is rather slack. Accordingly, when , can be guaranteed. But the parameter is usuallly a very small nonnegative number. This means that is a larger upper bound of the time-step.

(2) In this case, . (i)If and hold, then the numerical solution obeys the upper bound; that is, if and .(ii)If holds, then an arbitrary time-step can always insure if and .

*Proof. * is clear.

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 holds, the above is equivalent to
Thus, if and the conditions of (i) hold, the numerical solution according to Scheme (15) is always less than 1.

(ii) If and , the inequality always holds, and an arbitrary time-step can always insure .

*Remark 3. *In Scheme (15), if the decay rate of biomass is greater than the maximum specific growth rate , 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 is less than the recombination coefficient , then the biomass is growth that is governed by a nonlinear advection-diffusion equation.

#### 4. Numerical Examinations

We now provide some numerical examples to examine the performance of our method. Consider a nondimensional spatial domain . 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.

For comparison purposes, the following standard upwind forward Euler finite-difference method (EE) is also coded:

*Example 1. *Let and . We consider an initial mass density profile with the form of
where , , , , and ; , , , , and ; , , , , and . The following parameters are employed in our computations: , ; , , , and , . The homogeneous Neumann boundary conditions are imposed on all grid points on the boundary (4 edges) of . The results are plotted in Figure 1, using .

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 is larger than the maximum specific growth rate . 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 , 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 and . We consider an initial density profile with the form of
where , , and ; , , , , and ; , , , , and . The following parameters are employed in our computations: , ; , , , and , . The homogeneous Dirichlet boundary conditions are imposed on all grid points on the boundary (4 edges) of . The results are plotted in Figure 3, using .

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 is much smaller than the maximum specific growth rate . 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 , 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 and reach to 10^{−2} at 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 , but presented Scheme (15) can ensure that the numerical solution is nonnegative and bounded. Therefore, our Scheme (15) is stable for Example 2.

*Example 3. *Let and . We consider an initial density profile with the form of (34), where , , and ; , , and ; , , , , and . Fix the model parameters: , ; , , , and , . We use the same initial conditions as in Example 3. The results are plotted in Figure 4, using .

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 is smaller than the combined effective coefficient . Because the nutrition function is attenuating, which brings about the decrease of the recombination coefficient over time, the biomass decay rate will be larger than the combined effective coefficient . 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 , 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 and reach to 10^{−2} at 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 , but presented Scheme (15) can ensure that the numerical solutions are nonnegative and bounded. Therefore, our Scheme (15) is also stable for Example 3.

*Example 4. *Let and consider an initial profile of the biomass in the form of
The following parameters are employed in our computations: , and , using . Dirichlet conditions were imposed on the boundary of . The numerical results are plotted in Figure 9.

It can be observed from Figure 9 that the results are nonnegative and bounded for the problems with discontinuous initial value using presented Scheme (15).

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.

#### 5. 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 well-designed 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.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

The work by the 2nd author is partially supported by the United States NSF Grant under Award no. 1214188, the United States ARO Contract no. W911NF-12-1-0147, and the United States DoD (DTRA) Grant no. HDTRA1-13-1-0025. It is also partially supported by the Open Research Fund Program of the State Key Laboratory of Advanced Technology of Design and Manufacturing for Vehicle Body, Hunan University, China, under Grant no. 41215002. The work by the 3rd author is partially supported by the Natural Science Foundation for Young Scientists of Shanxi, China, Grant no. 2012021002-2.