Abstract
In this work, we introduce a linear finitedifference methodology to approximate nonnegative and bounded solutions of a coupled system of nonlinear parabolic partial differential equations that describes the growth of two different microbial colonies on a substrate of nutrients. Some simpler versions of the model of interest possess qualitative results that guarantee the existence and uniqueness of nonnegative and bounded solutions. However, the exact determination of analytical solutions of our system (and even of those simpler versions of our model) corresponding to physically meaningful initial conditions, may be a difficult task, whence the need to design computational methods to approximate the solutions of our mathematical model is pragmatically justified. Our numerical technique has the advantage of conditionally preserving the nonnegative and bounded characters of initialboundary data. The most important analytical results of this work are summarized as a theorem of existence and uniqueness of nonnegative and bounded numerical solutions, whose proof relies on the nonsingularity property of Mmatrices, and the fact that the entries of the inverses of these matrices are positive real numbers. We provide some illustrative simulations to evince the fact that our method preserves the mathematical characteristics of the solutions that we mentioned above.
MSC: 65Q10, 65M06, 39A14, 92D25.
Keywords:
computational modeling; biological films; positivity and boundedness; linear finitedifference scheme; food safety investigation1 Introduction
The problem of modeling mathematically the growth dynamics of biological films is an important topic of investigation in view of the many practical problems where these complex structures appear. For instance, microbial films find applications in the development of biological techniques to treat contaminated fluids [14], in the design of microbial fuel cells to produce electricity by means of biological or chemical procedures [58], in the production of new generations of sensors through bacterial signaling systems [9] and in the investigation of the corrosion of material surfaces for environmental engineering purposes [10,11], among many other scientific and engineering problems of pragmatic relevance [12,13]. Needless to mention that, according to a report by the Chinese National Institute of Health [14], biological films are likely to be directly or indirectly responsible for most of the bacterial infections in humans.
The mathematical modeling of biological films ideally offers the advantage of providing an inexpensive means to determine some analytical properties of the solutions of quantitative paradigms in general scenarios. However, even simple biological film structures possess a high level of complexity, which is difficult to be faithfully reflected in a single system of partial differential equations. As a consequence, most of the realistic, mathematical models used to describe the growth dynamics of microbial colonies are complex systems of equations for which the exact determination of meaningful analytical solutions is a difficult task. As expected, this complexity is substantially increased when one considers not only the interaction between a single colony of bacteria with the substrate, but also when different types of bacteria are present in the medium.
In the present work, we consider a nonlinear system of diffusive partial differential equations that models the interaction between two different types of microbial colonies and a substrate of nutrients. The mathematical model under investigation is a slightly more general version of the one proposed in [15], and it is an extension of the systems of differential equations studied in [1618]. Our model reflects many of the characteristics found experimentally in microbial films, like
(a) the existence of a sharp front of biomass at the fluidsolid transition,
(b) the presence of a threshold of biomass density,
(c) the fact that the spreading of biomass is significant only when the biomass is close to the threshold,
(d) the application of reaction kinetics mechanisms in the production of biomass,
(e) the compatibility of the biomass spreading mechanism with hydrodynamics and with nutrient transferconsumption models,
among other features observed in the practice. The functions of interest in our model (concentration of nutrients and population sizes) are nonnegative functions that have been normalized with respect to maximum allowed quantities, so that the sets of possible values that these variables may take on are subsets of .
In view of the difficulties to calculate analytical solutions of the model under study, we propose here a computational method to approximate them. Motivated by some previous and successful efforts in the design of numerical techniques to estimate the solutions of much simpler forms of the system of equations investigated in the present work [19,20], we follow a linear approach [2126] to provide a linear and implicit finitedifference discretization of the mathematical model of interest. After some algebraic manipulations, we will readily verify that our technique may be represented conveniently in vector form using a square matrix with real entries which, under suitable conditions, turns out to be an Mmatrix, that is, a strictly diagonally dominant matrix with nonpositive offdiagonal entries, and positive components in its main diagonal [27,28].
Following some approaches of the specialized literature employed for much simpler models [2932], the main properties of Mmatrices will be used in order to establish the conservation of the nonnegative and the bounded characters of approximations. More concretely, we will use the facts that every Mmatrix is nonsingular and that all the entries of their inverses are positive numbers [33]. Using these tools, we will prove a theorem of existence and uniqueness of nonnegative and bounded solutions of our technique. As expected, we will provide some simulations in order to evince the performance of our method and illustrate the facts that nonnegativity and boundedness are preserved at each iteration when some suitable conditions on the model and the computational parameters are satisfied.
The present work is divided as follows. In Section 2, we introduce the mathematical model under investigation, the physical relevance of the model parameters, the mathematical nomenclature and some conventions. Some particular models, for which theorems on the existence and uniqueness of nonnegative and bounded solutions are available, are described in that stage of our work. Section 3 introduces the numerical technique employed to approximate the solutions of our model. We introduce therein the discrete notation employed in this manuscript and provide a useful vector representation of our method. In Section 4, we derive some analytical properties of our finitedifference scheme (like its capability to conditionally preserve the nonnegative and bounded characters of approximations), and we present some illustrative simulations obtained through a computational implementation of our technique. Finally, we close this work with a section of concluding remarks.
2 Preliminaries
2.1 Mathematical model
Let p be a positive integer, and let Ω be a bounded subset of which is open and connected. Throughout this manuscript, will denote the closure of Ω in , while will represent the closure of the set of positive numbers in ℝ. In this work, s, u and v will be real functions defined in the set , which are twice differentiable in the interior of their domains and satisfy the coupled system of parabolic partial differential equations
for every . Here, the set Ω will represent a spatial domain which may physically be identified with a Petri dish for practical purposes, and the variable t denotes time. The spatial operators ∇, ∇⋅ and are, respectively, the gradient, the divergence and the Laplacian.
Physically, the system of equations (1) describes the dynamics of interaction between a substrate of nutrients s and two different biological colonies labeled here and , with respective relative densities u and v that have been normalized with respect to a maximum value of the biological mass. In system (1), the function w is the total biological mass of the system at each point of Ω and each instant of time, that is, for every ,
For the sake of concreteness, we may think of the growth dynamics of biological films formed of colonies of pseudomonas putida and listeria monocytogenes on a Petri dish with a substrate of nutrients. Moreover, in our investigation the diffusion factors are provided by the expressions
and the reaction functions are given by
In our model, the constants (the substrate diffusion coefficient), and (the diffusion coefficients corresponding to and , respectively), and (the maximum specific consumption rate of each colony), and (the Monad halfsaturation constant of each of the colonies), and (the maximum specific growth rate of anaerobic growth for and , respectively), and (the decay rate of the total biomass of the system for each colony), (the maximum specific growth rate of aerobic growth), α and β are all nonnegative numbers with .
On the other hand, J is a real matrix of size , which represents the direction of preferred biomass spreading across the substrate. In the twodimensional case, J is a 2by2 diagonal matrix of the form
where θ is a fixed number in the interval . In this case, one readily verifies that is the modified divergence operator
Before we close this section, it is worthwhile to notice that system (1) is a generalized form of the model investigated in [15]. In fact, our model reduces to the one studied in that work if we let . In the present manuscript, the distinction between the decay rate of the total biomass for each colony is made in order to derive a particular model of interaction between a substrate of nutrients and a single microbial colony.
2.2 Substratebiomass model
It is worthwhile to notice that (1) reduces to a previously studied model in the investigation of microbial films when the biomass is formed exclusively by one colony type. Indeed, letting v be equal to zero, the resulting system of equations may be rewritten as
where u is the biomass density normalized with respect to a maximum possible value, and s represents the density of the substrate of nutrients. Meanwhile, the function is reduced to
for every . Appropriate initialboundary conditions are required. For instance, one may consider
Let ℱ be the real function defined on through the expression
Our work is greatly motivated by the next analytical result.
Proposition 1Letandsatisfy the following conditions:
Then there exists a unique solution of system (11) subject to (13), which satisfies the following properties:
Proof The proof is a direct consequence of Theorems 2.1 and 2.2 of [18]. □
It is important to remark the fact that this result guarantees the existence and uniqueness of nonnegative and essentially bounded solutions of the particular biofilm model (11). However, even for this simplified version of (1), the calculation of particular solutions for nontrivial initial conditions is a complicated task. These observations motivate the numerical approach reported in this manuscript.
2.3 Simple model
Let be a continuous function. A further simplification of (1) leads to the nonlinear partial differential equation investigated in [19,20], which is a model that considers the presence of one colony of bacteria with no dynamical interaction with the substrate. The mathematical system investigated in those reports is described by the following equation:
for every . Obviously, suitable initial conditions must be imposed on Ω. Following the approach of Section 2.2, we consider initialboundary conditions of the form
Under these circumstances, the following result is an immediate corollary of Proposition 1. It establishes sufficient conditions for the existence and uniqueness of nonnegative and bounded solutions of Equation (15) subject to (16).
Corollary 1Letrandϕsatisfy:
Then there exists a unique solutionuof (15) subject to the initialboundary conditions (16), which satisfies the following properties:
3 Computational model
3.1 Nomenclature
For the rest of this work, we will restrict our attention to the case . Lower and higher dimensional scenarios are treated in a similar fashion.
Throughout this manuscript, we will assume that K, M and N are positive integers and fix the spatial domain of . We will take uniform partitions of , of the form
and
for every and . Clearly, the corresponding partition norms are given by
Fix also a temporal period of length equal to , and take a uniform partition of the interval of the form
be the corresponding partition norm. For each , and , we use the notation , and to represent approximations to the exact values of s, u and v, respectively, at the point . Moreover, define
In addition, we will employ the following standard forwarddifference operators, defined for , and for every , and :
Clearly, these linear operators provide firstorder approximations of the partial derivatives of q with respect to x, y and t, respectively, at the point . We also define the discrete operators
which are respectively approximations of order and of the secondorder partial derivative of s with respect to x and with respect to y. Let
We need to introduce now some nonstandard operators. Let , and let , and . To that end, define
Also, let
For any real number r, define
For the sake of convenience, we introduce the function given by the expression
3.2 Finitedifference scheme
Following the nomenclature introduced in Section 3.1, the finitedifference methodology employed here to approximate the solutions of the biological film model (1) is given by the system of recursive equations
for every , and . In addition, we will impose discrete, homogeneous Neumann boundary data for each of the three functions of interest.
The following remarks will be useful in Section 4.1 in order to establish the main properties of our method. They will also serve as an intermediate step in order to present our finitedifference scheme in vector form.
Remark 1 It important to notice that the three difference equations in (36) are approximations of the biological film model (1), which may be written in explicit form after we perform suitable algebraic steps. Indeed, for instance, the difference equation employed to approximate the dynamics of the substrate function may be rewritten as
where
for , and every , and . Obviously, the coefficients and are negative, the coefficient is positive, and the inequality holds.
We analyze now the second recursive formula of (36).
Remark 2 On the other hand, the discrete equation in (36), which approximates the dynamics of u, may be written as
where
for , , , and for every , and . Observe that the coefficients in Equation (41) are all nonpositive numbers and that the constants are all less than or equal to zero if and only if
holds. Observe that (45) is satisfied in the particular case when , and the approximation to the substrate function at the point and time is nonnegative. Finally, notice that the inequality
holds if and only if
Next, we examine the last iterative formula of (36).
Remark 3 Finally, we consider the last difference equation of (36), which approximates the values of the function v. In this case, it is also easy to verify that the finitedifference approximation is equivalent to the equation
for each , and . Here, the coefficients β are given by the expression (42), and
It is important to recall that the coefficients are nonnegative for , and that all the constant is nonpositive if and only if
is satisfied. This last inequality holds in particular when , and the approximation to the substrate function at the point and time is nonnegative. Finally, assuming that is nonpositive and is positive, then the inequality
holds if and only if
is satisfied.
3.3 Vector representation
The present section is devoted to present the new system of algebraic equations (38), (41) and (48) conveniently in a single and recursive vector representation. To that end, we arrange the approximations of the method at the time in lexicographic order. More precisely, we let
Let be the juxtaposition of the vectors , and in that order, and let be the lexicographically ordered vector of the approximate total biomass of the system at the time , that is, let
Throughout, we define the vector exactly as the vector , except that the components corresponding to boundary points are equal to zero.
Let k be an element of . We use the notation I to represent the identity matrix of size , and let be the diagonal matrix of size , all of whose entries in the diagonal are ones, except the first and the last. Define the tridiagonal matrix as the matrix of the same size of I given by
for every . Define the block matrix of size through
where the zeros are zero matrices of size .
On the other hand, for each , let be the matrix of the same size of given by
and let be the diagonal matrix of the same size as provided by
Define then the block matrix of the same size as via
Let be the diagonal matrix of size whose diagonal components are the constants arranged in lexicographic order for all those coordinates which do not correspond to the boundaries of Ω; otherwise, let those entries be equal to zero. Define the matrix in a similar fashion as , but using now the coefficients instead of .
Finally, for each , let be the matrix of the same size of provided by the formula
and introduce the block matrix of the same size as through
Remark 4 In view of the nomenclature introduced in this work and the algebraic presentation of the recursive equations of our method, the iterative system of equations (38), (41) and (48) may be presented in vector form as
where is the block matrix of size defined by
and is the lexicographically ordered vector of initial approximations. Clearly, the zeros of this matrix are zero matrices of size .
In view of these remarks, the method introduced in this work is linear and implicit, and provides linear approximations to the solutions of (1). Computationally, our technique is coded using an implementation of the stabilized biconjugate gradient method, which is a method that has been employed in several works to solve sparse systems arising in biology [19,20], among other disciplines of the natural science and engineering.
4 Results
In this section, we establish firstly some results on the existence and uniqueness of nonnegative and bounded solutions of the finitedifference method (36) using the vector representation (65). In a second stage, we provide some illustrative simulations that reflect the capability of the method to yield nonnegative and bounded approximations for nonnegative and bounded initial profiles.
4.1 Analytical results
The analytical cornerstone of our investigation is the concept of Mmatrices. Recall that an Mmatrix is a real, square matrix A which satisfies the following three properties:
(i) the offdiagonal elements of A are nonpositive numbers,
(ii) the diagonal entries are positive numbers, and
(iii) the matrix A is strictly diagonally dominant.
In general, the Mmatrices are important in many areas of mathematics; in particular, they are useful tools in the field of numerical analysis, where the fact that Mmatrices are nonsingular is employed extensively. Moreover, all the entries of the inverse of any Mmatrix are positive numbers [33]. These facts are summarized conveniently in the following result.
Lemma 4.1 (Fujimoto [33])
Every Mmatrix is invertible, and the entries of its inverse matrix are all positive numbers.
We say that a real vector x is nonnegative (respectively, positive) if all of its components are nonnegative (respectively, positive) numbers; this fact is represented by (respectively, ). We say that x is bounded from above (respectively, strictly bounded from above) by 1 if all the components of x are less than or equal to (respectively, strictly less than) 1, a fact that is denoted by (respectively, ). It is readily checked that holds if and only if , where e is the vector of the same size of x, all of whose components are equal to 1. We employ the notation (respectively, ) to signify that the vector x is nonnegative and bounded from above (respectively, strictly bounded from above) by 1.
The following result establishes that the matrix in Equation (66) is an Mmatrix under suitable conditions. We will use here the remarks of Section 3.2 and will follow the notation introduced previously. Throughout, k will be an element of the set .
Lemma 4.2Assume thatand. Thenis anMmatrix if all the conditions (46), (47), (51) and (53) hold for everyand.
Proof By Remarks 1, 2, 3 and the inequalities (46) and (51), the offdiagonal entries of are nonpositive. In turn, the inequalities (47) and (53) assure the strict diagonal dominance of this matrix. We conclude then that is an Mmatrix by definition. □
The following is the main analytical result of this work. It establishes sufficient conditions that guarantee the existence and uniqueness of nonnegative and bounded new approximations at each iteration of the finitedifference scheme (36).
Theorem 4.3 (Existence and uniqueness)
If, thenholds if the inequalities (46) and (51), as well as
and
are all satisfied for everyand.
Proof Beforehand, notice that the inequalities (67) and (68) imply respectively that (47) and (53) also hold. Under these circumstances, is an Mmatrix and, by Lemma 4.1, it is nonsingular and the entries of its inverse are all positive numbers. Also, and, thus, . It follows that . In order to establish the boundedness from above, let e be the vector of the same dimension as , all whose components are equal to 1, and let . Notice that Equation (65) becomes
where
It is important to notice that the rows of matrix associated to the homogeneous Neumann boundary conditions contain a 1 and a −1 as their only nonzero entries, and that the corresponding components of are equal to zero. In view of that, those entries of are equal to zero. We analyze next the rest of the rows of , focusing on each of the three block rows of the presentation (66).
• Row 1. Evidently, the corresponding component of is given by the expression for some and .
• Row 2. In this case, the corresponding component of the vector is
• Row 3. For each row of the third block row of the matrix , the corresponding component of assumes the form
We have verified that all the components of are nonnegative in all cases. Thus, or, equivalently, that is bounded from above by 1, as desired. □
4.2 Numerical examples
In this section, we provide some numerical simulations in order to evince the performance of our method in general and, in particular, its capability to preserve the nonnegative and bounded characters of solutions. Throughout, we will restrict our attention to the domain of . The boundary conditions will be homogeneous of the Neumann type, and the initial biomass profiles considered will take on the form
for every which belongs to Ω. Here, L is a suitable positive integer, , and are positive numbers, and are points in the interior of Ω for each . For the sake of convenience, we define the vectors for each such l.
Our first example provides a qualitative comparison against some of the results reported in [19,20]. To that end, observe that the partial differential equation (15) is obtained from system (1) by fixing the constants , , . In addition, we must consider the initial conditions and for every , and any fixed constant .
Example 1 Consider system (1) with the parameter values and initial conditions of the previous paragraph, in which case the model equation (15) results. Let the initial profile of the function u be given by (73), with ; , , , , ; , , , , ; , , , , . Fix the model constants , and , and the computational parameters and . Under these circumstances, Figures 14 show the approximate solution of the function u at the times 0, 6, 8 and 10, respectively. The results are in good qualitative agreement with respect to those reported in Figure 2 of [20], which correspond to Example 2.1 therein. Moreover, the present example illustrates the capability of our method the preserve the nonnegativity and the boundedness of the solutions.
Figure 1. Approximate solutionuat time0. Graph of the approximate solution u of Equation (1) on the domain at the time 0. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 1.
Figure 2. Approximate solutionuat time6. Graph of the approximate solution u of Equation (1) on the domain at the time 6. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 1.
Figure 3. Approximate solutionuat time8. Graph of the approximate solution u of Equation (1) on the domain at the time 8. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 1.
Figure 4. Approximate solutionuat time10. Graph of the approximate solution u of Equation (1) on the domain at the time 10. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 1.
For our following example, we study the dynamics of growth of a biological film whose interaction with the (nonconstant) substrate is described by (11). Observe that that system is obtained from (1) by letting , , and letting and being fixed positive constants, say, equal to 1. Also, we need to impose the initial condition for every .
Example 2 Fix the set of parameter values and initial data described in the preceding paragraph for the model (1) whence the simplified model (11) results. In addition, we consider the same initial profile for the function u employed in Example 1, and let for every . The rest of the parameters are , , , , , , and . Computationally, we used the values and . Under these conditions, Figures 58 show the approximate solution of u at the times 10, 25, 50 and 100, respectively. The insets are the corresponding approximations to the substrate function s.
Figure 5. Approximate solutionssanduat time2.5. Graphs of the approximate solutions s and u of Equation (1) on the domain at the time 12.5. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 2.
Figure 6. Approximate solutionssanduat time25. Graphs of the approximate solutions s and u of Equation (1) on the domain at the time 25. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 2.
Figure 7. Approximate solutionssanduat time50. Graphs of the approximate solutions s and u of Equation (1) on the domain at the time 50. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 2.
Figure 8. Approximate solutionssanduat time100. Graphs of the approximate solutions s and u of Equation (1) on the domain at the time 100. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 2.
We consider now the full system (1).
Example 3 Let , , , , , , , , , , , , and . Fix a homogeneous initial substrate equal to 0.9. The initial profile of the colony is given by (73) with ; , , , ; , , , ; , , , . Meanwhile, the colony has ; , , , ; , , , ; , , , . Computationally, we let , and fix . Figures 912 present the results of simulating the growth dynamics of the biofilm system (1) under these circumstances. The graphs show snapshots of the solutions of u and v at the times 2.5 and 25, while the insets provide the corresponding approximations of the substrate and the total biological mass of the system. The results give evidence of the complex dynamics of growth of the two colonies of bacteria and the substrate of nutrients. We have observed that the conditions of nonnegativity and boundedness are satisfied at each iteration, facts which are in agreement with Theorem 4.3.
Figure 9. Approximate solutionsuandsat time2.5. Graph of the approximate solution u of Equation (1) on the domain at the time 2.5. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 3. The inset is a checkerboard graph corresponding to the function s.
Figure 10. Approximate solutionsvandwat time2.5. Graph of the approximate solution v of Equation (1) on the domain at the time 2.5. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 3. The inset is a checkerboard graph corresponding to the function w.
Figure 11. Approximate solutionsuandsat time25. Graph of the approximate solution u of Equation (1) on the domain at the time 25. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 3. The inset is a checkerboard graph corresponding to the function s.
Figure 12. Approximate solutionsvandwat time25. Graph of the approximate solution v of Equation (1) on the domain at the time 25. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 3. The inset is a checkerboard graph corresponding to the function w.
Next, we approximate the solutions of the problem proposed in Example 3 at longer periods of time.
Example 4 As a continuation of the previous example, we consider the same model with the same parameters and initial condition. Figures 1316 present the approximate solutions at the times 50 and 100. Observe that the biomass has consumed most of the nutrients of the substrate by the time 50. After a period of time of length 100, the substrate function is basically identically equal to zero. It is worth noticing that the nonnegative character of the substrate and biomasses is preserved at each iteration, even during longer periods of time. These facts are in perfect agreement with Theorem 4.3. We have obtained simulations for longer periods of time. The results reflect that the solution u tends to reach a constant function as t goes to infinity, which is in obvious agreement with the expression of the second equation of system (1). Meanwhile, we have observed that v tends to increase approximately as a horizontal plane, as t increases. These observations are also in agreement with the third equation of our model.
Figure 13. Approximate solutionsuandsat time50. Graph of the approximate solution u of Equation (1) on the domain at the time 50. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 4. The inset is a checkerboard graph corresponding to the function s.
Figure 14. Approximate solutionsvandwat time50. Graph of the approximate solution v of Equation (1) on the domain at the time 50. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 4. The inset is a checkerboard graph corresponding to the function w.
Figure 15. Approximate solutionsuandsat time100. Graph of the approximate solution u of Equation (1) on the domain at the time 100. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 4. The inset is a checkerboard graph corresponding to the function s.
Figure 16. Approximate solutionsvandwat time100. Graph of the approximate solution v of Equation (1) on the domain at the time 100. The model and computational parameters, as well as the initial and boundary conditions, are given in Example 4. The inset is a checkerboard graph corresponding to the function w.
5 Conclusions
In this work, we designed a finitedifference method to approximate the solutions of a mathematical model that appears in the investigation of food safety. The model is a coupled system of parabolic partial differential equations with nonlinear diffusion and reaction laws, which describes the interaction between a substrate function and two different microbial colonies. For simpler versions of the problem, the specialized literature provides suitable theorems on the existence and uniqueness of solutions which are nonnegative and bounded from above at all times. However, even in the mathematically simpler versions of our system of equations, the exact determination of analytical solutions for significant initialvalue problems is a difficult task.
The method we introduced in this manuscript is a convenient, linear discretization of the mathematical model of interest. After some calculations, we showed that the method is a linear and implicit technique which may be rewritten in vector form through the multiplication by a real square matrix that, under suitable conditions on the parameters of the analytical and computational models, turns out to be an Mmatrix. The facts that every Mmatrix is nonsingular and that all the entries of their inverses are positive numbers are the most important tools to establish conditions that guarantee the existence and uniqueness of nonnegative and bounded approximations for every set of nonnegative and bounded initial conditions. The main result is summarized as a theorem of existence and uniqueness of numerical solutions, and the computer simulations provided show that the properties of nonnegativity and boundedness are preserved at every iteration of the method when the conditions of our result are satisfied.
Before we close this work, we must mention that we have obtained more computational simulations of the evolution of solutions of model (1) with various sets of values of the model parameters. On the one hand, the results evince the complex dynamics of this biological system, but they also motivate future numerical studies toward the resolution of practical problems in the investigation of biological films.
Competing interests
The author declares that he has no competing interests.
Authors’ contributions
The main author was the only contributor of this manuscript.
Acknowledgements
The author wants to express his deepest gratitude to Dr. I. E. MedinaRamírez of the Department of Chemistry at the Universidad Autónoma de Aguascalientes for all her invaluable comments and enlightening discussions during the preparation of this manuscript. Also, he wishes to thank the anonymous reviewers for their invaluable comments and criticisms.
References

Cao, W, Zhang, H, Wang, Y, Pan, JZ: Bioremediation of polluted surface water by using biofilms on filamentous bamboo. Ecol. Eng.. 42, 146–149 (2012)

Castillo, I, Hernández, P, Lafuente, A, RodríguezLlorente, ID, Caviedes, MA, Pajuelo, E: Selfbioremediation of corkprocessing wastewaters by (chloro) phenoldegrading bacteria immobilized onto residual cork particles. Water Res.. 46, 1723–1734 (2011). PubMed Abstract  Publisher Full Text

Ellwood, NTW, Di Pippo, F, Albertano, P: Phosphatase activities of cultured phototrophic biofilms. Water Res.. 46, 378–386 (2011). PubMed Abstract  Publisher Full Text

Jiao, Y, Zhao, Q, Hao, X, Jin, W, You, S: Bioaugmentation of a biological contact oxidation ditch with indigenous nitrifying bacteria for in situ remediation of nitrogenrich stream water. Bioresour. Technol.. 102, 990–995 (2011). PubMed Abstract  Publisher Full Text

Jain, A, Gazzola, G, Panzera, A, Zanoni, M, Marsili, E: Visible spectroelectrochemical characterization of geobacter sulfurreducens biofilms on optically transparent indium tin oxide electrode. Electrochim. Acta. 56, 10776–10785 (2011). Publisher Full Text

Liu, Y, Harnisch, F, Fricke, K, Schröder, U, Climent, V, Feliu, JM: The study of electrochemically active microbial biofilms on different carbonbased anode materials in microbial fuel cells. Biosens. Bioelectron.. 25, 2167–2171 (2010). PubMed Abstract  Publisher Full Text

Liu, Y, Harnisch, F, Fricke, K, Sietmann, R, Schröder, U: Improvement of the anodic bioelectrocatalytic activity of mixed culture biofilms by a simple consecutive electrochemical selection procedure. Biosens. Bioelectron.. 24, 1006–1011 (2008). PubMed Abstract  Publisher Full Text

Yang, Y, Sun, G, Guo, J, Xu, M: Differential biofilms characteristics of hewanella decolorationis microbial fuel cells under open and closed circuit conditions. Bioresour. Technol.. 102, 7093–7098 (2011). PubMed Abstract  Publisher Full Text

Checa, SK, Zurbriggen, MD, Soncini, FC: Bacterial signaling systems as platforms for rational design of new generations of biosensors. Curr. Opin. Biotechnol.. 23, 766–772 (2012). PubMed Abstract  Publisher Full Text

Little, BJ, Lee, JS, Ray, RI: The influence of marine biofilms on corrosion: a concise review. Electrochim. Acta. 54, 2–7 (2008). Publisher Full Text

Upadhyayula, VKK, Gadhamshetty, V: Appreciating the role of carbon nanotube composites in preventing biofouling and promoting biofilms on material surfaces in environmental engineering: a review. Biotechnol. Adv.. 28, 802–816 (2010). PubMed Abstract  Publisher Full Text

Kim, MJ, Baek, S, Jung, SH, Cho, KH: Dynamical characteristics of bacteria clustering by selfgenerated attractants. Comput. Biol. Chem.. 31, 328–334 (2007). PubMed Abstract  Publisher Full Text

Velusamy, V, Arshak, K, Korostynska, O, Oliwa, K, Adley, C: An overview of foodborne pathogen detection: in the perspective of biosensors. Biotechnol. Adv.. 28, 232–254 (2010). PubMed Abstract  Publisher Full Text

Chen, HH, Liu, X, Ni, C, Lu, YP, Xiong, GY, Lu, YY, Wang, SQ: Bacterial biofilms in chronic rhinosinusitis and their relationship with inflammation severity. Auris, Nasus, Larynx. 39, 169–174 (2012). PubMed Abstract  Publisher Full Text

Eberl, HJ, Schraft, H: A diffusionreaction model of a mixedculture biofilm arising in food safety studies. Mathematical Modeling of Biological Systems, pp. 109–120. Springer, Berlin (2008)

Eberl, HJ, Demaret, L: A finite difference scheme for a degenerated diffusion equation arising in microbial ecology. Electron. J. Differ. Equ.. 15, 77–95 (2007)

Eberl, HJ, Parker, D, van Loosdrecht, M: A new deterministic spatiotemporal continuum model for biofilm development. J. Theor. Med.. 3, 161–176 (2001). Publisher Full Text

Efendiev, MA, Eberl, HJ, Zelik, SV: Existence and longtime behavior of solutions of a nonlinear reactiondiffusion system arising in the modeling of biofilms. RIMS Kokyuroko. 1258, 49–71 (2002)

MoralesHernández, MD, MedinaRamírez, IE, AvelarGonzález, FJ, MacíasDíaz, JE: An efficient recursive algorithm in the computational simulation of the bounded growth of biological films. Int. J. Comput. Methods. 9, (2012) Article ID 1250050

MoralesHernández, MD, MedinaRamírez, IE, AvelarGonzález, FJ, MacíasDíaz, JE: Corrigendum: “An efficient recursive algorithm in the computational simulation of the bounded growth of biological films”. Int. J. Comput. Methods. 10, (2013) Article ID 1392001

Mickens, RE, Jordan, PM: A positivitypreserving nonstandard finite difference scheme for the damped wave equation. Numer. Methods Partial Differ. Equ.. 20, 639–649 (2004). Publisher Full Text

Mickens, RE, Jordan, PM: A new positivitypreserving nonstandard finite difference scheme for the DWE. Numer. Methods Partial Differ. Equ.. 21, 976–985 (2005). Publisher Full Text

Mickens, RE: Nonstandard finitedifference schemes for reactiondiffusion equations. Numer. Methods Partial Differ. Equ.. 15, 201–214 (1999). Publisher Full Text

Mickens, RE: Discretizations of nonlinear differential equations using explicit nonstandard methods. J. Comput. Appl. Math.. 110, 181–185 (1999). Publisher Full Text

Ongun, MY, Arslan, D, Garrappa, R: Nonstandard finite difference schemes for fractional order Brusselator system. Adv. Differ. Equ.. 2013, (2013) Article ID 102

Liao, C, Ding, X: Nonstandard finite difference variational integrators for nonlinear Schrödinger equation with variable coefficients. Adv. Differ. Equ.. 2013, (2013) Article ID 12

Ding, D, Ma, Q, Ding, X: A nonstandard finite difference scheme for an epidemic model with vaccination. J. Differ. Equ. Appl.. 19, 179–190 (2011)

Zhang, Z, Chen, D: A new criterion on existence and uniqueness of stationary distribution for diffusion processes. Adv. Differ. Equ.. 2013, (2013) Article ID 13

MacíasDíaz, JE: A bounded finitedifference discretization of a twodimensional diffusion equation with logistic nonlinear reaction. Int. J. Mod. Phys. C. 22, 953–966 (2011). Publisher Full Text

MacíasDíaz, JE, Puri, A: A boundednesspreserving finitedifference scheme for a damped nonlinear wave equation. Appl. Numer. Math.. 60, 934–948 (2010). Publisher Full Text

MacíasDíaz, JE, RuizRamírez, J: A nonstandard symmetrypreserving method to compute bounded solutions of a generalized NewellWhiteheadSegel equation. Appl. Numer. Math.. 61, 630–640 (2011). Publisher Full Text

RuizRamírez, J, MacíasDíaz, JE: A finitedifference scheme to approximate nonnegative and bounded solutions of a FitzHughNagumo equation. Int. J. Comput. Math.. 88, 3186–3201 (2011). Publisher Full Text

Fujimoto, T, Ranade, RR: Two characterizations of inversepositive matrices: The HawkinsSimon condition and the Le ChatelierBraun principle. Electron. J. Linear Algebra. 11, 59–65 (2004)