Mathematical Problems in Engineering
Volume 2015 (2015), Article ID 708497, 17 pages
http://dx.doi.org/10.1155/2015/708497

Research Article

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

G. F. Sun,¹ G. R. Liu,1,2 and M. Li¹

¹College of Mathematics, Taiyuan University of Technology, China
²School of Aerospace Systems, University of Cincinnati, 2851 Woodside Dr, Cincinnati, OH 45221, USA

Received 16 July 2014; Accepted 29 October 2014

Academic Editor: Seungik Baek

Copyright © 2015 G. F. Sun et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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⁻³ for Example 1.(2) is smaller and macro error is smaller. This means that our model is convergent for variable .

Table 1

Spatial mesh

Time-step

Macro error
CPU time (s)

0.02
0.00004
0.1
0.0017
0.0015
0.0014

0.02
0.00005
0.125
0.0019
0.0017
0.0016

0.02
0.0001
0.25
0.0030
0.0028
0.0027
857.3280

0.02
0.00016
0.4
0.0043
0.0041
0.0040
536.7030

0.02
0.0002
0.5
0.0052
0.0050
0.0049
425.5460

0.02
0.00032
0.8
0.0078
0.0076
0.0075
259.5790

0.02
0.0004
1.0
0.0095
0.0094
0.0093
209.6880

0.02
0.0005
1.25
0.0117
0.0116
0.0115
169.1560

Figure 1: Numerical solutions of Example 1 governed by (1) on at four different times, , and 10, using finite-difference Scheme (15) with . The initial profile is given by (33).

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

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⁻³ at and reach to 10⁻² 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.

Table 2

Spatial mesh

Time-step

Macro error
CPU time (s)

0.02
0.00004
0.1

0.0164

0.02
0.00005
0.125

0.0181

0.02
0.0001
0.25

0.0016
0.0237

0.02
0.00016
0.4

0.0025
0.0280
563.2970

0.02
0.0002
0.5

0.0031
0.0302
435.2030

0.02
0.00032
0.8

0.0045
0.0349
263.7350

0.02
0.0004
1.0

0.0054
0.0371
224.6560

0.02
0.0005
1.25

0.0064
0.0394
172.6880

Figure 3: Numerical solutions of Example 2 governed by equations set (1) on at four different times, , and , using finite-difference Scheme (15) with . The initial profile is given by (34).

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

Figure 5: The numerical results of Scheme (15) and EE method, using the same initial-boundary condition and calculated parameters as 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⁻³ at and reach to 10⁻² 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.

Table 3

Spatial mesh

Time-step

Macro error
CPU time (s)

0.02
0.00004
0.1

0.0438

0.02
0.00008
0.2

0.0486

0.02
0.00016
0.4
0.001
0.0540

692.3440

0.02
0.0002
0.5
0.0012
0.0558
0.0010
595.8440

0.02
0.00032
0.8
0.0018
0.0596
0.0020
315.2340

0.02
0.0004
1.0
0.0021
0.0618
0.0027
295.6090

Figure 6: Numerical solutions of Example 3 governed by equations set (1) on at four different times, , and 12, using finite-difference Scheme (15) with . The initial profile is given by (34).

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

Figure 8: Numerical solutions of Example 3 governed by equations set (1) on at four different times, , and 12, using the EE method with . The initial profile is given by (34).

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

Figure 9: Numerical solutions of Example 4 governed by equations set (1) on at four different times, , and 10, using finite-difference Scheme (15) with . The initial profile is given by (35).

Figure 10: The comparison figures of Scheme (15) and EE method. Numerical solutions of Example 4 governed by equations set (1) on at two different times, , 2.4. The initial profile is given by (35). The model parameters: ; , , , and ; (, ).

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.

References

H. J. Eberl, D. F. Parker, and M. C. M. Vanloosdrecht, “A new deterministic spatio-temporal continuum model for biofilm development,” Journal of Theoretical Medicine, vol. 3, no. 3, pp. 161–175, 2001.

M. A. Efendiev, H. J. Eberl, and S. V. Zelik, “Existence and longtime behavior of solutions of a nonlinear reaction-diffusion system arising in the modeling of biofilms,” RIMS Kokyuroko, no. 1258, pp. 49–71, 2002.

T. Fujimoto and R. R. Ranade, “Two characterizations of inverse-positive matrices: the Hawkins-Simon condition and the Le Chatelier-Braun principle,” Electronic Journal of Linear Algebra, vol. 11, pp. 59–65, 2004.

H. A. Abdusalam, “Multiple soliton solutions for the Nagumo equation and the modified general Burgers-Fisher equation,” International Journal of Computational Methods, vol. 3, no. 3, pp. 371–381, 2006.

H. A. Abdusalam and E. S. Fahmy, “Traveling wave solutions for nonlinear wave equation with dissipation and nonlinear transport term through factorizations,” International Journal of Computational Methods, vol. 4, no. 4, pp. 645–652, 2007.

H. J. Eberl and L. Demaret, “A finite difference scheme for a degenerated diffusion equation arising in microbial ecology,” Electronic Journal of Differential Equations, vol. 15, pp. 77–95, 2007.

C. R. Armbruster, T. S. Forster, R. M. Donlan, H. A. O'Connell, A. M. Shams, and M. M. Williams, “A biofilm model developed to investigate survival and disinfection of Mycobacterium mucogenicum in potable water,” Biofouling, vol. 28, no. 10, pp. 1129–1139, 2012.

M. R. Frederick, C. Kuttler, B. A. Hense, and H. J. Eberl, “A mathematical model of quorum sensing regulated EPS production in biofilm communities,” Theoretical Biology and Medical Modelling, vol. 8, no. 1, article 8, 2011.

Q. Liao, Y.-J. Wang, Y.-Z. Wang et al., “Two-dimension mathematical modeling of photosynthetic bacterial biofilm growth and formation,” International Journal of Hydrogen Energy, vol. 37, no. 20, pp. 15607–15615, 2012.

Y.-Z. Wang, Q. Liao, X. Zhu, R. Chen, C.-L. Guo, and J. Zhou, “Bioconversion characteristics of Rhodopseudomonas palustris CQK 01 entrapped in a photobioreactor for hydrogen production,” Bioresource Technology, vol. 135, pp. 331–338, 2013.

Y. Imamura, J. Chandra, P. K. Mukherjee et al., “Fusarium and Candida albicans biofilms on soft contact lenses: model development, influence of lens type, and susceptibility to lens care solutions,” Antimicrobial Agents and Chemotherapy, vol. 52, no. 1, pp. 171–182, 2008.

J. E. Macías-Díaza, S. Macíasa, and I. E. Medina-Ramírezb, “An efficient nonlinear finite-difference approach in the computational modeling of the dynamics of a nonlinear diffusion-reaction equation in microbial ecology,” Computational Biology and Chemistry, vol. 47, pp. 24–30, 2013.

M. D. Morales-Hern{\'a}ndez, I. E. Medina-Ram{\'{\i}}rez, F. J. Avelar-Gonz{\'a}lez, and J. E. Mac{\'{\i}}as-D{\'{\i}}az, “An efficient recursive algorithm in the computational simulation of the bounded growth of biological films,” International Journal of Computational Methods, vol. 9, no. 4, 1250050, 15 pages, 2012.

M. R. Parsek and P. K. Singh, “Bacterial biofilms: an emerging link to disease pathogenesis,” Annual Review of Microbiology, vol. 57, pp. 677–701, 2003.

A. D. Peacock, Y.-J. Chang, J. D. Istok et al., “Utilization of microbial biofilms as monitors of bioremediation,” Microbial Ecology, vol. 47, no. 3, pp. 284–292, 2004.

R. Chen, Y.-K. Pu, Q. Liao, X. Zhu, and Y.-Z. Wang, “A simulation on PSB biofilm formation with considering cell inactivation,” International Journal of Hydrogen Energy, vol. 38, no. 35, pp. 15670–15679, 2013.

G. Roeselers, M. C. M. V. Loosdrecht, and G. Muyzer, “Phototrophic biofilms and their potential applications,” Journal of Applied Phycology, vol. 20, no. 3, pp. 227–235, 2008.

J. Ruiz-Ramíreza and J. E. Macías-Díaza, “A finite-difference scheme to approximate non-negative and bounded solutions of a FitzHugh-Nagumo equation,” International Journal of Computer Mathematics, vol. 88, no. 15, pp. 3186–3201, 2011.