Research

# On a linear finite-difference model of a mixed-culture biological system arising in food safety studies

Jorge E Macías-Díaz

Author Affiliations

Departamento de Matemáticas y Física, Universidad Autónoma de Aguascalientes, Avenida Universidad 940, Cd. Universitaria, Aguascalientes, 20131, Mexico

Advances in Difference Equations 2013, 2013:357  doi:10.1186/1687-1847-2013-357

 Received: 22 July 2013 Accepted: 18 November 2013 Published: 5 December 2013

This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### Abstract

In this work, we introduce a linear finite-difference methodology to approximate non-negative 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 non-negative 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 non-negative and bounded characters of initial-boundary data. The most important analytical results of this work are summarized as a theorem of existence and uniqueness of non-negative and bounded numerical solutions, whose proof relies on the non-singularity property of M-matrices, 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 finite-difference scheme; food safety investigation

### 1 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 [1-4], in the design of microbial fuel cells to produce electricity by means of biological or chemical procedures [5-8], 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 [16-18]. Our model reflects many of the characteristics found experimentally in microbial films, like

(a) the existence of a sharp front of biomass at the fluid-solid 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 transfer-consumption models,

among other features observed in the practice. The functions of interest in our model (concentration of nutrients and population sizes) are non-negative 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 [21-26] to provide a linear and implicit finite-difference 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 M-matrix, that is, a strictly diagonally dominant matrix with non-positive off-diagonal entries, and positive components in its main diagonal [27,28].

Following some approaches of the specialized literature employed for much simpler models [29-32], the main properties of M-matrices will be used in order to establish the conservation of the non-negative and the bounded characters of approximations. More concretely, we will use the facts that every M-matrix is non-singular 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 non-negative 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 non-negativity 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 non-negative 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 finite-difference scheme (like its capability to conditionally preserve the non-negative 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

(1)

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 ,

(2)

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

(3)

(4)

and the reaction functions are given by

(5)

(6)

(7)

(8)

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 half-saturation 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 non-negative 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 two-dimensional case, J is a 2-by-2 diagonal matrix of the form

(9)

where θ is a fixed number in the interval . In this case, one readily verifies that is the modified divergence operator

(10)

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 Substrate-biomass 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

(11)

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

(12)

for every . Appropriate initial-boundary conditions are required. For instance, one may consider

(13)

for suitable functions .

Let ℱ be the real function defined on through the expression

(14)

Our work is greatly motivated by the next analytical result.

Proposition 1Letandsatisfy the following conditions:

(A) andfor every,

(B) and,

(C) for every, and.

Then there exists a unique solution of system (11) subject to (13), which satisfies the following properties:

1. ,

2. ,

3. for every, and.

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 non-negative and essentially bounded solutions of the particular biofilm model (11). However, even for this simplified version of (1), the calculation of particular solutions for non-trivial 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:

(15)

for every . Obviously, suitable initial conditions must be imposed on Ω. Following the approach of Section 2.2, we consider initial-boundary conditions of the form

(16)

Under these circumstances, the following result is an immediate corollary of Proposition 1. It establishes sufficient conditions for the existence and uniqueness of non-negative and bounded solutions of Equation (15) subject to (16).

Corollary 1Letrandϕsatisfy:

(A) andfor every,

(B) and,

(C) for every, and.

Then there exists a unique solutionuof (15) subject to the initial-boundary conditions (16), which satisfies the following properties:

1. ,

2. ,

3. for every, and.

### 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

(17)

and

(18)

for every and . Clearly, the corresponding partition norms are given by

(19)

(20)

Fix also a temporal period of length equal to , and take a uniform partition of the interval of the form

(21)

for every . Let

(22)

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

(23)

In addition, we will employ the following standard forward-difference operators, defined for , and for every , and :

(24)

(25)

(26)

Clearly, these linear operators provide first-order approximations of the partial derivatives of q with respect to x, y and t, respectively, at the point . We also define the discrete operators

(27)

(28)

which are respectively approximations of order and of the second-order partial derivative of s with respect to x and with respect to y. Let

(29)

We need to introduce now some non-standard operators. Let , and let , and . To that end, define

(30)

Also, let

(31)

(32)

where but . Let

(33)

For any real number r, define

(34)

For the sake of convenience, we introduce the function given by the expression

(35)

#### 3.2 Finite-difference scheme

Following the nomenclature introduced in Section 3.1, the finite-difference methodology employed here to approximate the solutions of the biological film model (1) is given by the system of recursive equations

(36)

for every , and . In addition, we will impose discrete, homogeneous Neumann boundary data for each of the three functions of interest.

Let , and let

(37)

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 finite-difference 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

(38)

where

(39)

(40)

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

(41)

where

(42)

(43)

(44)

for , , , and for every , and . Observe that the coefficients in Equation (41) are all non-positive numbers and that the constants are all less than or equal to zero if and only if

(45)

holds. Observe that (45) is satisfied in the particular case when , and the approximation to the substrate function at the point and time is non-negative. Finally, notice that the inequality

(46)

holds if and only if

(47)

whenever and .

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 finite-difference approximation is equivalent to the equation

(48)

for each , and . Here, the coefficients β are given by the expression (42), and

(49)

(50)

It is important to recall that the coefficients are non-negative for , and that all the constant is non-positive if and only if

(51)

is satisfied. This last inequality holds in particular when , and the approximation to the substrate function at the point and time is non-negative. Finally, assuming that is non-positive and is positive, then the inequality

(52)

holds if and only if

(53)

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

(54)

(55)

(56)

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

(57)

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

(58)

for every . Define the block matrix of size through

(59)

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

(60)

and let be the diagonal matrix of the same size as provided by

(61)

Define then the block matrix of the same size as via

(62)

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

(63)

and introduce the block matrix of the same size as through

(64)

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

(65)

where is the block matrix of size defined by

(66)

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 bi-conjugate 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 non-negative and bounded solutions of the finite-difference 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 non-negative and bounded approximations for non-negative and bounded initial profiles.

#### 4.1 Analytical results

The analytical cornerstone of our investigation is the concept of M-matrices. Recall that an M-matrix is a real, square matrix A which satisfies the following three properties:

(i) the off-diagonal elements of A are non-positive numbers,

(ii) the diagonal entries are positive numbers, and

(iii) the matrix A is strictly diagonally dominant.

In general, the M-matrices are important in many areas of mathematics; in particular, they are useful tools in the field of numerical analysis, where the fact that M-matrices are non-singular is employed extensively. Moreover, all the entries of the inverse of any M-matrix are positive numbers [33]. These facts are summarized conveniently in the following result.

Lemma 4.1 (Fujimoto [33])

Every M-matrix is invertible, and the entries of its inverse matrix are all positive numbers.

We say that a real vector x is non-negative (respectively, positive) if all of its components are non-negative (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 non-negative and bounded from above (respectively, strictly bounded from above) by 1.

The following result establishes that the matrix in Equation (66) is an M-matrix 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 anM-matrix 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 off-diagonal entries of are non-positive. In turn, the inequalities (47) and (53) assure the strict diagonal dominance of this matrix. We conclude then that is an M-matrix by definition. □

The following is the main analytical result of this work. It establishes sufficient conditions that guarantee the existence and uniqueness of non-negative and bounded new approximations at each iteration of the finite-difference scheme (36).

Theorem 4.3 (Existence and uniqueness)

If, thenholds if the inequalities (46) and (51), as well as

(67)

and

(68)

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 M-matrix and, by Lemma 4.1, it is non-singular 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

(69)

where

(70)

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

(71)

for suitable and .

• Row 3. For each row of the third block row of the matrix , the corresponding component of assumes the form

(72)

for some and .

We have verified that all the components of are non-negative 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 non-negative 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

(73)

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 1-4 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 non-negativity 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 (non-constant) 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 5-8 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 9-12 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 non-negativity 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 13-16 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 non-negative 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 finite-difference 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 non-negative 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 initial-value 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 M-matrix. The facts that every M-matrix is non-singular 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 non-negative and bounded approximations for every set of non-negative 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 non-negativity 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. Medina-Ramí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

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

2. Castillo, I, Hernández, P, Lafuente, A, Rodríguez-Llorente, ID, Caviedes, MA, Pajuelo, E: Self-bioremediation of cork-processing wastewaters by (chloro) phenol-degrading bacteria immobilized onto residual cork particles. Water Res.. 46, 1723–1734 (2011). PubMed Abstract | Publisher Full Text

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

4. 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 nitrogen-rich stream water. Bioresour. Technol.. 102, 990–995 (2011). PubMed Abstract | Publisher Full Text

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

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

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

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

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

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

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

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

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

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

15. Eberl, HJ, Schraft, H: A diffusion-reaction model of a mixed-culture biofilm arising in food safety studies. Mathematical Modeling of Biological Systems, pp. 109–120. Springer, Berlin (2008)

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

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

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

19. Morales-Hernández, MD, Medina-Ramírez, IE, Avelar-González, FJ, Macías-Dí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

20. Morales-Hernández, MD, Medina-Ramírez, IE, Avelar-González, FJ, Macías-Dí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

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

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

23. Mickens, RE: Nonstandard finite-difference schemes for reaction-diffusion equations. Numer. Methods Partial Differ. Equ.. 15, 201–214 (1999). Publisher Full Text

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

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

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

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

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

29. Macías-Díaz, JE: A bounded finite-difference discretization of a two-dimensional diffusion equation with logistic nonlinear reaction. Int. J. Mod. Phys. C. 22, 953–966 (2011). Publisher Full Text

30. Macías-Díaz, JE, Puri, A: A boundedness-preserving finite-difference scheme for a damped nonlinear wave equation. Appl. Numer. Math.. 60, 934–948 (2010). Publisher Full Text

31. Macías-Díaz, JE, Ruiz-Ramírez, J: A non-standard symmetry-preserving method to compute bounded solutions of a generalized Newell-Whitehead-Segel equation. Appl. Numer. Math.. 61, 630–640 (2011). Publisher Full Text

32. Ruiz-Ramírez, J, Macías-Díaz, JE: A finite-difference scheme to approximate non-negative and bounded solutions of a FitzHugh-Nagumo equation. Int. J. Comput. Math.. 88, 3186–3201 (2011). Publisher Full Text

33. Fujimoto, T, Ranade, RR: Two characterizations of inverse-positive matrices: The Hawkins-Simon condition and the Le Chatelier-Braun principle. Electron. J. Linear Algebra. 11, 59–65 (2004)