- Email: [email protected]
Analytical substrate flux approximation for the Monod boundary value problem
Applied Mathematics and Computation 218 (2011) 1484–1494
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Analytical substrate flux approximation for the Monod boundary value problem Fazal Abbas, Hermann J. Eberl ⇑ Dept. Mathematics and Statistics, University of Guelph, Guelph, ON, Canada
a r t i c l e
i n f o
Keywords: Two point boundary value problem Monod kinetics Biofilm model Analytical approximation
a b s t r a c t We present an analytical approximation for the diffusive flux of a substrate into a reactive layer, in which the substrate is degraded according to Monod kinetics. This problem is described by a nonlinear two-point boundary value problem. The approximation is derived based on a Homotopy Analysis Method idea and verified computationally, by comparison against a numerical solution of the problem. The analytical approximation is easy to evaluate and depends only on model parameters. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction We are concerned with the computation of the diffusive flux of a dissolved substrate into a homogeneous reactive layer, where the substrate is degraded according to Monod kinetics. The motivation for our study is the mathematical modeling for bacterial biofilms, but problems of this type arise in other applications in Mathematical Biology as well, e.g. in tumor modeling, etc. Bacterial biofilms are microbial depositions on submerged interfaces. Bacteria adhere to the surface (called substratum in the biofilm literature) and start to produce extracellular polymeric substances (EPS), in which they are themselves embedded, and which gives them protection against harmful environmental impacts, such as mechanical washout, or exposure to biocides. Microbial growth in biofilms requires the availability of sufficient nutrients and, in the aerobic case, oxygen. These dissolved substrates enter the biofilm from the surrounding aqueous phase. In the biofilm, they diffuse and are degraded by the microbes. Due to the disparity of the characteristic time-scales of bacterial growth and substrate decay, both processes are often decoupled and the transport-reaction problem for the dissolved substrate can be studied in a quasi-steady state, assuming a biofilm of fixed thickness [19]. In traditional one-dimensional biofilm models, this leads to a semi-linear two point boundary value problem. Due to the steady state assumption, microbial growth in the biofilm can be equated with the substrate consumption in the biofilm, which can be equated with the flux of substrate into the biofilm. For given parameters, the substrate concentration in the biofilm and thus the flux from the bulk liquid can be computed conveniently with numerical methods. However, from a modeling and analysis point of view, it would often be desirable to have a simple closed form algebraic expression for the substrate flux, that depends only on model parameters. This plays an important role, for example, in upscaling of the model as well as for analytical treatment, e.g. the study of long term behavior of biofilm growth models, because it would allow to reduce the model from a free boundary value problem to an ordinary differential equation. In this study we will present a simple algebraic approximation of the substrate flux into a single-species, single-substrate biofilm. This approximation is derived using arguments from Homotopy Analysis/Perturbation Methods. It is verified by comparison against a numerical solution of the two-point boundary value problem. ⇑ Corresponding author. E-mail address: [email protected] (H.J. Eberl). 0096-3003/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.05.102
1485
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
Perturbation methods are an extensively used technique to derive approximate solutions for nonlinear differential equations. Their downside is the assumption that small parameters must exist, which sometimes limits the method’s applicability. To compensate for this limitation various other approximation methods have been developed such as Homotopy Analysis Method (HAM), Homotopy Perturbation Method (HPM) and Adomian’s decomposition methods. The idea of homotopy methods stems from Differential Topology [9]. The underlying approach of HAM and HPM is that the solution of the nonlinear problem can be expressed as a series of solutions of simpler problems and that this series can be truncated after a finite number of terms to give a sufficiently accurate approximation. In HAM, this series is a generalized Taylor series in a homotopy parameter p 2 [0, 1] which satisfies the original (not elementary solvable) problem exactly for p = 1 and a simpler, solvable problem for p = 0. In HPM this is a perturbation series in p. Homotopy Perturbation Methods have been extensively used to solve two point-boundary value problems and they work generally well for linear problems or problems in which the nonlinearity is polynomial or allows for an unconditionally convergent series expansion, e.g. in [3–7,14,18]. However, it has been pointed out that in other cases convergence problems can arise [1,8]. It was suggested that the Homotopy Analysis Method can be understood as a generalization of the Homotopy Perturbation Method that often can overcome these difficulties, at the expense of introducing an additional slack variable, commonly referred to as convergence control parameter [11,12]. However, for the flux calculation problem at hand we can find an algebraically simpler, efficient way to deal with such convergence problems. 2. Problem formulation The traditional one-dimensional steady state biofilm model was proposed in [15], where the growth of biomass is controlled by the substrate flux across the water-biofilm interface. It is assumed that the biomass is homogeneously distributed in the biofilm and that the dissolved substrate follows the two-point boundary value problem with Monod reaction kinetics for the dependent variable substrate concentration c [M L3]
0¼D
@ 2 c lX 1 @z2 Y
c
jþc
@c ð0Þ ¼ 0; @z
;
cðkÞ ¼ cb ;
ð1Þ
where l is the maximum specific growth rate [T1], D is the diffusion coefficient [L2 T1], X1 is the biomass density [M L3], j is the Monod Half Saturation Concentration [M L3] and Y is a yield coefficient [–]. One set of parameters from the biofilm literature [16] is presented in Table 1. These are the values that we shall use in our calculations. We denote by k the biofilm thickness [L], and by cb the bulk substrate concentration [M L3]. The boundary condition at z = 0 models that the substratum is impermeable to substrate, while the boundary condition at z = k prescribes the substrate concentration at the water/biofilm interface. The solution of (1) can be obtained analytically if j cb or if j cb as a solution to limiting linear problems. In the general case, the analytical solution is not known. From an application point of view, we are in particular interested in the diffusive substrate flux Ddc into the biofilm. dz z¼k Integrating (1) once gives
Jðcb ; kÞ :¼
dc lX 1 ðkÞ ¼ dz YD
Z
k
0
c dz; jþc
ð2Þ
where we choose our notation as to indicate that the flux depends on both, the bulk substrate concentration cb and the biofilm thickness k. We non-dimensionalize (1) using the scalings x ¼ kz ; s ¼ jc and sb ¼ cjb and obtain 2
d s 2
dx
¼ ak2
where a ¼
X 1l . YDj
s ; 1þs
sð1Þ ¼ sb ;
ds ð0Þ ¼ 0; dx
ð3Þ
Instead of (2) we obtain
jðsb ; kÞ :¼
Z 1 ds s 2 dx: ¼ ak dxx¼1 1þs 0
ð4Þ
Our aim is to find a closed form analytical approximation of (4).
Table 1 Model parameters from [16]. Symbol
l D X1
j Y
Parameter Maximum growth rate Diffusion coefficient Biomass density Monod constant Yield coefficient
Units 1
T L2 T1 M L3 M L3 –
Value 1.65 105 s1 1.0 105 cm 2 s1 2.5 103 g cm3 3.9 106 g cm3 0.071
1486
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
3. Analytical series solution 3.1. Derivation Starting point is (3), re-written in the form 2
2
d s
d s þ s 2 ak2 s ¼ 0; 2 dx dx
ds ð0Þ ¼ 0: dx
sð1Þ ¼ sb ;
ð5Þ
We define for p 2 [0, 1] the homotopy
"
v ðx; pÞ
H1 ðx; pÞ ¼ ð1 pÞ
# " # @2 @2 @2 @2 2 v ðx; pÞ v ðx; pÞ v ðx; 0Þ þ p v ðx; pÞ þ v ðx; pÞ v ðx; pÞ ak v ðx; pÞ ¼ 0; @x2 @x2 @x2 @x2
d v ð0; pÞ ¼ 0: dx
v ð1; pÞ ¼ sb ;
ð6Þ
It is easy to verify that for p = 1 this is equivalent with (3). Moreover, every function that satisfies the boundary conditions is a valid solution of (6) for p = 0. We expand v(x; p) in a Taylor series about p = 0 to obtain
v ðx; pÞ ¼ v ðx; 0Þ þ
þ1 X
v ðmÞ ðx; pÞjp¼0 m!
m¼0
!
pm :
ð7Þ
Formally, provided there is sufficient regularity in the problem and the series is convergent, we can write, therefore,
v ðx; 1Þ ¼ v ðx; 0Þ þ
þ1 X
v ðmÞ ðx; pÞjp¼0
m¼0
m!
ð8Þ
:
We call the terms under the sum the mth order deformation derivative, defined as
sm ðxÞ :¼
v ðmÞ ðx; pÞjp¼0 m!
ð9Þ
:
The solution of (5) is then
sðxÞ ¼ s0 þ
þ1 X
sm ðxÞ;
ð10Þ
m¼1
where s0(x) = v(x; 0). The material flux into the biofilm can be computed from this according to the first equality in (4) by differentiation and multiplication with a constant. The deformation derivatives can be computed in analogy to the theory summarized in [10]. For our specific problem we have the following algorithm. Theorem 3.1. For s0(x) = sb, the deformation derivatives sm associated with homotopy H1 are obtained recursively as solutions of the boundary value problems m X k¼0
2
sk
d smk 2
dx
2
þ
d sm1 2
dx
ak2 sm1 ¼ 0;
sm ð1Þ ¼ 0;
dsm ð0Þ ¼ 0: dx
ð11Þ
Proof. Substituting (7) into (6) and applying Leibniz’ Rule, we obtain for the left hand side of (6)
m X m L:H:S ¼ ð1 pÞðlÞ l l¼0
"
@2 @2 v ðx; pÞ 2 v ðx; pÞ v ðx; pÞ 2 v ðx; 0Þ @x @x
#ðmlÞ
or
L:H:S ¼
m X m k¼0
k
v ðkÞ ðx; pÞ
m1 Xm 1 @ 2 ðmkÞ @ 2 ðmk1Þ ðkÞ v ðx; pÞj m v ðx; pÞ v ðx; pÞjp¼0 : p¼0 @x2 @x2 k k¼0
For the right hand side we obtain
" #ðmlÞ m 2 X m @2 ðlÞ @ 2 R:H:S ¼ v ðx; pÞ þ v ðx; pÞ 2 v ðx; pÞ ak v ðx; pÞ ðpÞ @x2 @x l l¼0
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
1487
or
m1 X m1 @ 2 ðm1Þ @2 v ðx; pÞj m v ðkÞ ðx; pÞ 2 v ðm1kÞ ðx; pÞjp¼0 þ mak2 v ðm1Þ ðx; pÞjp¼0 : p¼0 2 @x @x k k¼0
R:H:S ¼ m
Rewriting this in terms of the deformation derivatives sm(x) this becomes
R:H:S ¼ m!
" 2 d sm1 2
dx
þ
m1 X
#
2
sk
d sm1k
ak2 sm1
2
dx
k¼0
and
" L:H:S ¼ m!
m X
2
d smk
sk
2
dx
k¼0
m1 X
2
d smk1
sk
2
dx
k¼0
# :
Thus, for the mth deformation derivative we obtain with the boundary conditions in Eq. (6) the linear boundary value problem m X
2
sk
k¼0
2
d smk dx
þ
2
d sm1 2
dx
ak2 sm1 ¼ 0;
dsm ð0Þ ¼ 0: dx
sm ð1Þ ¼ 0;
Note that for given s1 ; . . . ; sm1 , the function sm can be computed by integrating twice. Since we choose s0(x) = sb = const, the mth order derivative sm is an even polynomial of degree less than or equal 2m. Thus, s(x) is a conventional power series for the solutions of (3). Its coefficients can be computed conveniently with computer algebra software. We use Maple for this task. We now consider the homotopy
"
# " # @2 @2 @2 @2 2 H2 ðx; pÞ ¼ ð1 pÞ v ðx; pÞ 2 v ðx; 0Þ þ p 2 v ðx; pÞ þ v ðx; pÞ 2 v ðx; pÞ ak v ðx; pÞ ¼ 0; @x2 @x @x @x d v ð0; pÞ ¼ 0: dx
v ð1; pÞ ¼ 1;
ð12Þ
Theorem 3.2. For s0(x) = sb, the deformation derivatives sm associated with homotopy H1 are obtained recursively as solutions of the boundary value problems 2
d sm 2
dx
¼
m1 X
2
sk
k¼0
d
s 2 m1k
dx
þ ak2 sm1 ;
sm ð1Þ ¼ 0;
dsm ð0Þ ¼ 0: dx
Proof. As in the previous proof, we substitute (7) into (12) and use Leibniz’s rule to obtain
" #ðmlÞ m 2 X m @2 ðlÞ @ L:H:S ¼ v ðx; pÞ 2 v ðx; 0Þ ð1 pÞ @x2 @x l l¼0 or
L:H:S ¼
@ 2 ðmÞ @ 2 ðm1Þ v ðx; pÞj m v ðx; pÞjp¼0 : p¼0 @x2 @x2
Moreover, we obtain
" #ðmlÞ m 2 X m @ 2 v ðx; pÞ ðlÞ @ v ðx; pÞ 2 R:H:S ¼ þ v ðx; pÞ ak v ðx; pÞ ðpÞ @x2 @x2 l l¼0 or
R:H:S ¼ m
m1 X m1 @ 2 ðm1Þ @ 2 ðm1kÞ ðkÞ v ðx; pÞj m v ðx; pÞ v ðx; pÞjp¼0 þ mak2 v ðm1Þ ðx; pÞjp¼0 : p¼0 @x2 @x2 k k¼0
Re-writing this in terms of the deformation derivative, we obtain
" R:H:S ¼ m!
2
d
s 2 m1
dx
þ
m1 X k¼0
#
2
sk
d
s 2 m1k
dx
2
ak sm1
ð13Þ
1488
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
and
" L:H:S ¼ m!
# 2 d : s s 2 m 2 m1 dx dx 2
d
Thus, 2
d sm 2
dx
¼
m1 X k¼0
sk
d
2
dx
s 2 m1k
þ ak2 sm1 ;
sm ð1Þ ¼ 0;
dsm ð0Þ ¼ 0: dx
Again one obtains sm(x) as an even polynomial of degree less than or equal 2m. 3.2. Numerical comparison The recursions in Theorems 3.1 and 3.2 allow us to compute approximations to the solutions of (3) but do not give any insight how fast, if at all, these approximations converge. Every finite order approximation is a polynomial approximation of the solution. One problem that may arise is that the power series representation of the Monod nonlinearity in (3) has only a finite radius of convergence which indicates a finite radius of convergence for the polynomial approximation of the solution as well. This was accounted for to some extent by providing two homotopies, one to be used for sb < 1 and one for sb > 1. By construction, the coefficients of sm(x) according to Theorems 3.1 and 3.2 depend on the parameters ak2 and sb. Therefore, it P can be expected that these affect the convergence of the series sðxÞ ¼ 1 m¼0 sm ðxÞ. Due to the algebraic complexity, a rigorous convergence proof is difficult. Moreover, mere convergence of the series does not automatically imply that truncating the series after few terms leads to a good quantitative approximation of the solution of the underlying two-point boundary value problem. Therefore, we carry out a computational comparison of finite-term approximations according to (3.1) and (3.2) with a numerical solution. The latter is computed with Maple’s built in Finite Difference method. In fact, since our original goal was to determine the flux (4) of substrate into the reactive layer, we do not compare s(x) but the approximate flux s0 (1) against the numerically computed one. The set of model parameters in this study is taken from [16] and listed in Table 1. The lth order fluxes according to the analytical approximation are computed for k and cb = sbj as
J l ðcb ; kÞ ¼
l jX
k
s0m ð1Þ;
ð14Þ
m¼0
where s0m ð1Þ is calculated with the computer algebra software Maple according to the recursion in Theorem 3.1 for cb < j and according to Theorem 3.2 for cb > j. The results are plotted in Fig. 1 for four selected values of cb and for selected flux approximations Jl. In all cases we observe an acceptable agreement of Jl with the numerical solution J num ¼ dc ðkÞ for small dz k<~ k 2 105 cm only. For larger values of k, this approximation breaks down, with the possible exception of the largest bulk concentration value tested. From an application point of view the small range of k for which the homotopy based methods work renders the approximation unusable. Bacterial cells have a diameter of order of k = 104 cm, i.e. biofilms, as bacterial layers, are considerably thicker than ~ k. Although this is a negative result, we can use the non- or slow-convergent homotopy based approximations to motivate an easy to evaluate flux approximation in the next section. 4. Quadratic approximation to the series solutions Upon inspection of the deformation derivatives sm(x) for m 6 13 according to Theorem 3.1 and for i 6 10 according to Theorem 3.2, we propose a new approximation. It is based on the observation that sm(x) can be expressed as a quadratic polynomial in x plus Oðk4 Þ terms. We treat the cases of Theorem 3.1 (used for sb > 1) and Theorem 3.2 (used for sb < 1) separately but will see that both will lead to the same compact formula for a flux approximation. This will be compared against the numerical solution. 4.1. Derivation Lemma 4.1. The solution sm(x) of the generalized mth boundary value problem (11) for m P 1 has the form
sm ðxÞ ¼
ð1Þm1 2 2 ak x 1 þ Oðk4 Þ: 2sm1 b
Proof. We show the assertion for m = 1 explicitly and use mathematical induction for the general case. The boundary value problem (11) for m = 1 is
ð15Þ
1489
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494 8e-09
2e-10
6e-09
J [g/cm44]
1.5e-10
J [g/cm4]
numerical J2 J=9 J13
numerical J2 J=8 J10
1e-10
4e-09
2e-09
5e-11
0
0 0
1e-05
2e-05
3e-05
4e-05
0
5e-05
1e-05
lambda [cm]
6e-09
8e-10
2e-05
3e-05
4e-05
3e-05
4e-05
lambda [cm]
numerical J2 J=9 J13
numerical J2 J=8 J10
6e-10
J [g/cm4]
J [g/cm4]
4e-09
4e-10
2e-09 2e-10
0
0 0
1e-05
2e-05
3e-05
4e-05
5e-05
0
lambda [cm]
1e-05
2e-05
lambda [cm]
Fig. 1. Comparison of flux approximation according to (14) and Theorem 3.2 for sb = cb/j < 1 (left column) and Theorem 3.1 sb = cb/j > 1 (right column). Note: j = 3.9 106 g/cm3, see Table 1.
s0 s001 þ s1 s000 þ s000 ak2 s0 ¼ 0;
s1 ð1Þ ¼ 0;
s01 ð0Þ ¼ 0:
Since s0(x) = sb, simple integration gives
s1 ðxÞ ¼
ak2 2 ðx 1Þ; 2
which is of form (15). Assume that the assertion holds for l = 1, . . . , m 1. From (15) we obtain
s00m1 ðxÞ ¼
ð1Þm2 2 ak þ Oðk4 Þ: sm2 b
Substituting sm1 ðxÞ and s00m1 ðxÞ into (11) leads to the integration problem
s0 s00m ðxÞ þ s00m1 ðxÞ þ Oðk4 Þ ¼ 0;
sm ð1Þ ¼ 0;
s0m ð0Þ ¼ 0:
Thus
sm ðxÞ ¼
ð1Þm1 2 2 ak ðx 1Þ þ Oðk4 Þ: 2sm1 b
Proposition 4.2. The analytical approximate solution s(x) of the Monod boundary value problem (3) according to homotopy (6) for sb > 1 has the form
1490
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
sðxÞ ¼ sb þ
ak2 ðx2 1Þ sb þ Oðk4 Þ: 2 1 þ sb
ð16Þ
Moreover,
s0 ð1Þ ¼
sb ak2 þ Oðk4 Þ: 1 þ sb
ð17Þ
Proof. With the previous lemma we obtain
sðxÞ ¼ sb þ
1 ak2 ðx2 1Þ X ð1Þm1 þ Oðk4 Þ: 2 sm1 b m¼0
For sb > 1, this is
sðxÞ ¼ sb þ
ak2 x2 1 sb þ Oðk4 Þ: 2 1 þ sb
which is in the required form. For the derivative at k = 1 it follows
s0 ð1Þ ¼
sb ak2 þ Oðk4 Þ: 1 þ sb
Thus, by dropping the Oðk4 Þ terms in Eq. (16) a candidate for a closed form analytical solution is obtained for the bulk concentration sb > 1. We repeat the same analysis for the solution s(x) obtained with homotopy (12). Lemma 4.3. The solution sm(x) of the generalized mth boundary value problem (13) for m P 1 has the form
sm ðxÞ ¼
ak2 ð1 x2 Þ 4 ð1Þm sm b þ Oðk Þ: 2
ð18Þ
Proof Again we use mathematical induction. For m = 1 we obtain
s0 s000 þ s001 ak2 s0 þ Oðk4 Þ ¼ 0;
s1 ð1Þ ¼ sb ;
s01 ð0Þ ¼ 0
and thus
s1 ¼
ak2 ðx2 1Þ sb þ Oðk4 Þ: 2
This is the same expression obtained in the previous case. The assertion is assumed to hold for l = 1, 2, . . . , m 1, which means
sm1 ðxÞ ¼
ak2 ð1 x2 Þ þ Oðk4 Þ ð1Þm1 sm1 b 2
is true. The second order derivative is then
s00m1 ðxÞ ¼
ak2 þ Oðk4 Þ: ð1Þm sm1 b 2
Substituting sm1 and s00m1 into (13) we obtain
s00m ¼ s0 s00m1 ;
sm ð1Þ ¼ 0;
sm ð0Þ ¼ 0
from which
sm ðxÞ ¼
ak2 ð1 x2 Þ 4 ð1Þm sm b þ Oðk Þ: 2
Proposition 4.4. The analytical approximate solution s(x) of the Monod boundary value problem (3) according to homotopy (12) for sb < 1 has the form
sðxÞ ¼ sb þ
ak2 ðx2 1Þ sb þ Oðk4 Þ: 2 1 þ sb
ð19Þ
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
1491
Moreover,
s0 ð1Þ ¼
sb ak2 þ Oðk4 Þ: 1 þ sb
ð20Þ
Proof. Using sm(x) from (18) gives
sðxÞ ¼ sb þ
1 ak2 ð1 x2 Þ X 4 ð1Þm sm b þ Oðk Þ: 2 m¼1
Thus, for sb < 1,
sðxÞ ¼ sb þ
ak2 ðx2 1Þ sb þ Oðk4 Þ: 2 1 þ sb
The derivative at k = 1 is then obtained as
s0 ð1Þ ¼
sb ak2 þ Oðk4 Þ: 1 þ sb
Again, dropping the Oðk4 Þ terms in (19) we derive a candidate for a closed form approximation of the solution of the Monod boundary value problem. We point out that the approximation candidates obtained for Theorem 3.1 (for sb > 1) and Theorem 3.2 (for sb < 1) take the same expression, which is continuous at sb = 1. 4.2. Numerical verification We repeat the numerical verification for the flux approximation obtained in Propositions 4.2 and 4.4. After dropping Oðk4 Þ term and converting back into physical units, we obtain for the derivative of the substrate concentration at the boundary
J analytical ¼
lX 1 YD
k
cb
j þ cb
ð21Þ
:
We compare this approximation against the numerical solution of (1), which we calculated with a standard finite difference method on a very fine grid. Converting the analytical expression back into physical units, we obtain
J num ¼
dc ðkÞ; dz
where the derivative is approximated by a finite difference from the concentration values in the last two grid points. We use the parameters from [16], summarized in Table 1. The biofilm thickness is varied over the range 20 lm < k < 200 lm. The bulk substrate concentration ranges over 0 6 cb 6 2 104 g/cm3. These are realistic ranges for many biofilm applications. For three selected biofilm thicknesses, k = 74 lm, 110 lm, 164 lm, the results are plotted in Fig. 2. The analytical approximation is found to be in good agreement with the numerical solution, in particular for smaller k and large cb. The quality of the approximation decreases as k becomes large for small cb. Qualitatively the same results were obtained with a second set of parameters, taken from [19] (data not shown). In all cases the analytical approximation is slightly larger than the numerical solution. This is in accordance with the following theoretical result. Proposition 4.5. The analytical approximation of the flux according to (17) or (20) after dropping the Oðk4 Þ terms overestimates the flux of the exact solution of (3). Proof. We first show that the Oðk4 Þ approximation ~s according to (16),
~sðxÞ ¼ sb þ
ak2 ðx2 1Þ sb 2 1 þ sb
is a lower estimate for s(x), i.e. that ~sðxÞ 6 sðxÞ. Indeed, we have on the one hand
~s00 ¼ ak2
sb 1 þ sb
and on the other hand
ak2
1 1 þ 12 ak2 ðx2 1Þ 1þs ~sðxÞ sb b : ¼ ak2 s 2 1 1 þ ~sðxÞ 1 þ sb 1 þ 2 ak ðx2 1Þ b 2 ð1þsb Þ
1492
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
0.001
0.0008
J [g/cm4]
0.0006
0.0004
0.0002 lambda=0.0164cm lambda=0.0110cm lambda=0.0074cm
0 0
5e-05
0.0001
0.00015
0.0002
Cb [g/cm3] Fig. 2. Analytical vs numerical computation of the specific substrate flux J for three values of biofilm thickness k as a function of the bulk concentration cb: The analytical solution is drawn as solid line, the numerical solution is represented by marks.
Since
sb 1þsb
< 1, the numerator in the last expression is larger than the denominator, thus
~s00 6 ak2
~sðxÞ : 1 þ ~sðxÞ
Hence, by standard comparison theorems for two-point boundary value problems, e.g. as in [2], we have
~sðxÞ 6 sðxÞ: With ~sð1Þ ¼ sð1Þ ¼ sb , it follows that
~s0 ð1Þ P s0 ð1Þ:
5. Simple applications In order to demonstrate how an analytical flux approximation can simplify biofilm models we include two examples. 5.1. Stability of biofilm growth The standard single-species/single substrate biofilm model [19] can be written in the form
dk ¼ dt
Z 0
k
lX 1
c kl X 1 dz dðkÞk; jþc
ð22Þ
where kl is the cell death rate and d(k) the rate of biofilm detachment. The substrate concentration c is given by (1). Together (1) and (22) form a free boundary value problem for c. From an application point of view, however, usually one is more interested either in the biofilm thickness k or in the amount of substrate converted by the biofilm, i.e. the flux of substrate into the biofilm. Using (21) and sb ¼ cjb , the problem is simplified to the ordinary differential equation
dk cb k ðkl þ dðkÞk ¼ FðkÞ; ¼A dt j þ cb
ð23Þ
where the constant A contains the original model parameters and does not depend on cb and k. It represents the growth rate of the biofilm and accounts for diffusion of the substrate in the biofilm. 0 b We compute F 0 ðkÞ ¼ jAc kl dðkÞ d ðkÞk. The stability of the model, therefore, depends on the parameters and the speþcb cific detachment rate function. The most commonly used detachment function is the so-called L2 detachment model, in which case d(k) = dk [19]. In this case the biofilm model (22), (1) becomes with our flux approximation the logistic equation
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
dk ¼ kðB dkÞ; dt
B :¼
Acb
j þ cb
1493
kl :
If the compounded net growth rate B is positive, a biofilm will establish, i.e. the equilibrium k = 0 is unstable and the equilibrium k = B/d is asymptotically stable. We note that such a logistic growth model was used ad hoc to be fitted against experimental data for detachment experiments in [17]. Our results give an interpretation of the coefficient B that allow us to relate it to the bulk substrate concentration and hence to environmental conditions. It is well understood and was demonstrated in experiments [13,17] that in fact the detachment rate depends on the reactor hydrodynamics. More specifically, it can be assumed that the detachment rate d(k) depends on the shear rate, and thus on the reactor type. A frequently used laboratory reactor is the so-called roto-torque reactor, where biofilm grows between two concentric cylinders, one of which rotates. This can be approximated by biofilm growing between two parallel plates, a distance D apart. One of the plates is stationary, the other one is assumed to move at a constant velocity u. As long as u is small enough that no Couette–Taylor instabilities occur, the flow profile is linear with slope u/(D k). Thus, the detachment rate is in this case proportional to u/(D k) [20]. Hence, from (22), (1) for the roto torque reactor with shear rate dependent detachment rate we have
dk c ¼k B ¼: FðkÞ; dt Dk where c is proportional to the reactor’s rotation speed. We have
F 0 ðkÞ ¼ B
ck c : D k ðD kÞ2
The trivial equilibrium k = 0 (no biofilm, washout) is stable, i.e. F0 (k) < 0 if the flow rate, i.e. c, is fast enough, or if the growth rate, i.e. B, is small enough. This is the case if cb is small. Depending on parameters, a second, positive, equilibrium k⁄ exists if 0 < k :¼ D Bc, i.e. if the trivial equilibrium is unstable. Here we can assume that B > 0 because otherwise no biofilm growth would be possible, even in the absence of detachment. If it exists, k⁄ is asymptotically stable. 5.2. Reactor scale analysis Similar as [13], we consider a vessel of volume V that is filled with biofilm carrier chips that provide additional surface area for biofilms to form. This mimics a Moving Bed Biofilm Reactor of the type used in wastewater engineering. The inflow and outflow rate to and from the vessel is Q. On inflow, the substrate concentration is c1. In the vessel, we assume that the substrate is completely mixed. The system is modeled by (22) and (1), together with a reactor mass balance for the substrate concentration in the bulk cb
dcb Q dc ¼ ðc1 cb Þ E ðkÞ; V dz dt where the new parameter E subsumes reaction–diffusion and reactor parameters. The last term models the consumption of substrate by the biofilm. Using (21) as before, the model reduces to the two-dimensional ordinary differential equation
dk cb k ¼A ðkl þ dðkÞÞk; dt j þ cb dcb cb k ; ¼ aðc1 cb Þ b dt j þ cb
ð24Þ ð25Þ
where parameters a, b subsume the other model parameters. This model can be studied in the phase plane with elementary techniques after the detachment rate function d(k) is specified. For example, for the standard choice d(k) = dk as above, the Jacobian of this system is obtained as
2 F ðk; cb Þ ¼ 4 0
Acb
jþcb
kl 2dk jbcþcb
b
3
Ajk ðjþcb Þ2
j a ðjbk þc Þ2
5:
b
c1 The system permits a trivial equilibrium k ; cb ¼ ð0; c1 Þ, for which we obtain the eigenvalues k1 ¼ A jþc c and 1 k2 ¼ a < 0. Thus, this equilibrium is unstable, i.e. a biofilm will be able to establish itself, if the bulk substrate concentration c1 c1 is high enough that the maximum growth rate is bigger than the cell death rate, i.e. A jþc > c, otherwise it is stable. In the 1 ⁄⁄ ⁄⁄ former case, a positive non-trivial equilibrium (k , c ) is obtained as the unique intersection of the strictly increasing funccb tion ^ k ¼ Ad jþc cd with the strictly decreasing function ~ k ¼ ab ðC 1 cb Þðj þ cb Þ. Then k⁄⁄ > 0 and c b < c1 . The Jacobian in this b case becomes 0
F k
; c b
2 ¼4
dk
3
Ajk ðjþc Þ2 b
j
bk Ab ðc þ dk Þ a ðjþc Þ2 b
5
1494
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
with negative trace and positive determinant. Thus k ; c is asymptotically stable. b Hence, in all three applications that we briefly sketched here, the analytical flux approximation allows us to reduce a hybrid ODE/PDE model to a simple ordinary differential equation that can be studied with elementary techniques. 6. Conclusion Two different homotopy equations (6) and (12) were proposed to obtain series solutions for the Monod two point boundary value problem. Both lead to approximations that are irrelevant for practical purposes. However, using that the biofilm thickness k as a small parameter enters these approximations to the 4th power, we were able to obtain an easy to evaluate Oðk4 Þ analytical approximation of the infinite order approximations, which was shown to be an accurate approximation of the numerically computed substrate flux into the biofilm. This approximation will allow to simplify modeling tasks, such as upscaling of biofilm models and stability analysis, by reducing model complexity. Acknowledgement This research was supported in parts by the Higher Education Commission of Pakistan through a doctoral scholarship for FA. References [1] S. Abbasbandy, The application of the homotopy analysis method to nonlinear equations arising in heat transfer, Phys. Lett. A 360 (2006) 109–113. [2] S. Carl, S. Heikkilä, Nonlinear Differential Equations in Ordered Spaces, Chapman & Hall, 2000. [3] M.S.H. Chowdhury, I. Hashim, Analytical solutions to heat transfer equations by homotopy-perturbation method revisited, Phys. Lett. A 372 (2008) 1240–1243. [4] M. Duman, Asymptotic expansions for the Sturm–Liouville problem by homotopy perturbation method, Appl. Math. Comput. 216 (2010) 492–496. [5] D.D. Ganji, The application of He’s homotopy perturbation method to nonlinear equations arising in heat transfer, Phys. Lett. A 355 (2006) 337–341. [6] J.H. He, Comparison of homotopy perturbation method and homotopy analysis method, Appl. Math. Comput. 156 (6) (2004) 527–539. [7] J.H. He, Homotopy perturbation method for solving boundary value problems, Phys. Lett. A 350 (2006) 87–88. [8] S. Liang, D.J. Jeffrey, Comparison of homotopy analysis method and homotopy perturbation method through an evolution equation, Commun. Nonlinear Sci. Numer. Simul. 14 (12) (2009) 4057–4064. [9] S.J. Liao, Beyond Perturbation: Introduction to The Homotopy Analysis Method, Chapman and Hall, 2004. [10] S.J. Liao, Notes on the homotopy analysis method: some definitions and theorems, Commun. Nonlinear Sci. Numer. Simul. 14 (4) (2009) 983–997. [11] S.J. Liao, An optimal homotopy-analysis approach for strongly nonlinear differential equations, Commun. Nonlinear Sci. Numer. Simul. 15 (8) (2010) 2003–2016. [12] S.J. Liao, On the relationship between the homotopy analysis method and Euler transfer, Commun. Nonlinear Sci. Numer. Simul. 15 (6) (2010) 1421– 1431. [13] A. Masic, J. Bengtsson, M. Christensson, Measuring and modeling the oxygen profile in a nitrifying Moving Bed Biofilm Reactor, Math. Biosci. 227 (1) (2010) 1–11. [14] A. Rafiq, M. Ahmed, S.A. Hussain, A general approach to specific second order ordinary differential equations using homotopy perturbation method, Phys. Lett. A 372 (2008) 4973–4976. [15] B.E. Rittmann, P.L. McCarty, Model of steady state biofilm kinetics, Biotechnol. Bioeng. 22 (11) (1980) 2343–2357. [16] C.T. Skowlund, Effects of biofilm growth on a steady state biofilm model, Biotechnol. Bioeng. 35 (5) (1990) 502–510. [17] Y.P. Tsai, Impact of flow velocity on the dynamic behaviour of biofilm bacteria, Biofouling 21 (5) (2005) 267–277. [18] Y.G. Wang, H.F. Song, D. Li, Solving two-point boundary value problems using combined homotopy perturbation method and Green’s function method, Appl. Math. Comput. 212 (2009) 366–376. [19] O. Wanner, H. Eberl, E. Morgenroth, D. Noguera, C. Picioreanu, B. Rittmann, M.C.M. van Loosdrecht, Mathematical Modeling of Biofilms, IWA Publishing, 2006. [20] T. Xu, Numerical and statistical analysis of detachment forces in biofilm colonies, M.Sc. Thesis, Dept. Mathematics and Statistics, University of Guelph, 2005.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Analytical substrate flux approximation for the Monod boundary value problem Fazal Abbas, Hermann J. Eberl ⇑ Dept. Mathematics and Statistics, University of Guelph, Guelph, ON, Canada
a r t i c l e
i n f o
Keywords: Two point boundary value problem Monod kinetics Biofilm model Analytical approximation
a b s t r a c t We present an analytical approximation for the diffusive flux of a substrate into a reactive layer, in which the substrate is degraded according to Monod kinetics. This problem is described by a nonlinear two-point boundary value problem. The approximation is derived based on a Homotopy Analysis Method idea and verified computationally, by comparison against a numerical solution of the problem. The analytical approximation is easy to evaluate and depends only on model parameters. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction We are concerned with the computation of the diffusive flux of a dissolved substrate into a homogeneous reactive layer, where the substrate is degraded according to Monod kinetics. The motivation for our study is the mathematical modeling for bacterial biofilms, but problems of this type arise in other applications in Mathematical Biology as well, e.g. in tumor modeling, etc. Bacterial biofilms are microbial depositions on submerged interfaces. Bacteria adhere to the surface (called substratum in the biofilm literature) and start to produce extracellular polymeric substances (EPS), in which they are themselves embedded, and which gives them protection against harmful environmental impacts, such as mechanical washout, or exposure to biocides. Microbial growth in biofilms requires the availability of sufficient nutrients and, in the aerobic case, oxygen. These dissolved substrates enter the biofilm from the surrounding aqueous phase. In the biofilm, they diffuse and are degraded by the microbes. Due to the disparity of the characteristic time-scales of bacterial growth and substrate decay, both processes are often decoupled and the transport-reaction problem for the dissolved substrate can be studied in a quasi-steady state, assuming a biofilm of fixed thickness [19]. In traditional one-dimensional biofilm models, this leads to a semi-linear two point boundary value problem. Due to the steady state assumption, microbial growth in the biofilm can be equated with the substrate consumption in the biofilm, which can be equated with the flux of substrate into the biofilm. For given parameters, the substrate concentration in the biofilm and thus the flux from the bulk liquid can be computed conveniently with numerical methods. However, from a modeling and analysis point of view, it would often be desirable to have a simple closed form algebraic expression for the substrate flux, that depends only on model parameters. This plays an important role, for example, in upscaling of the model as well as for analytical treatment, e.g. the study of long term behavior of biofilm growth models, because it would allow to reduce the model from a free boundary value problem to an ordinary differential equation. In this study we will present a simple algebraic approximation of the substrate flux into a single-species, single-substrate biofilm. This approximation is derived using arguments from Homotopy Analysis/Perturbation Methods. It is verified by comparison against a numerical solution of the two-point boundary value problem. ⇑ Corresponding author. E-mail address: [email protected] (H.J. Eberl). 0096-3003/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.05.102
1485
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
Perturbation methods are an extensively used technique to derive approximate solutions for nonlinear differential equations. Their downside is the assumption that small parameters must exist, which sometimes limits the method’s applicability. To compensate for this limitation various other approximation methods have been developed such as Homotopy Analysis Method (HAM), Homotopy Perturbation Method (HPM) and Adomian’s decomposition methods. The idea of homotopy methods stems from Differential Topology [9]. The underlying approach of HAM and HPM is that the solution of the nonlinear problem can be expressed as a series of solutions of simpler problems and that this series can be truncated after a finite number of terms to give a sufficiently accurate approximation. In HAM, this series is a generalized Taylor series in a homotopy parameter p 2 [0, 1] which satisfies the original (not elementary solvable) problem exactly for p = 1 and a simpler, solvable problem for p = 0. In HPM this is a perturbation series in p. Homotopy Perturbation Methods have been extensively used to solve two point-boundary value problems and they work generally well for linear problems or problems in which the nonlinearity is polynomial or allows for an unconditionally convergent series expansion, e.g. in [3–7,14,18]. However, it has been pointed out that in other cases convergence problems can arise [1,8]. It was suggested that the Homotopy Analysis Method can be understood as a generalization of the Homotopy Perturbation Method that often can overcome these difficulties, at the expense of introducing an additional slack variable, commonly referred to as convergence control parameter [11,12]. However, for the flux calculation problem at hand we can find an algebraically simpler, efficient way to deal with such convergence problems. 2. Problem formulation The traditional one-dimensional steady state biofilm model was proposed in [15], where the growth of biomass is controlled by the substrate flux across the water-biofilm interface. It is assumed that the biomass is homogeneously distributed in the biofilm and that the dissolved substrate follows the two-point boundary value problem with Monod reaction kinetics for the dependent variable substrate concentration c [M L3]
0¼D
@ 2 c lX 1 @z2 Y
c
jþc
@c ð0Þ ¼ 0; @z
;
cðkÞ ¼ cb ;
ð1Þ
where l is the maximum specific growth rate [T1], D is the diffusion coefficient [L2 T1], X1 is the biomass density [M L3], j is the Monod Half Saturation Concentration [M L3] and Y is a yield coefficient [–]. One set of parameters from the biofilm literature [16] is presented in Table 1. These are the values that we shall use in our calculations. We denote by k the biofilm thickness [L], and by cb the bulk substrate concentration [M L3]. The boundary condition at z = 0 models that the substratum is impermeable to substrate, while the boundary condition at z = k prescribes the substrate concentration at the water/biofilm interface. The solution of (1) can be obtained analytically if j cb or if j cb as a solution to limiting linear problems. In the general case, the analytical solution is not known. From an application point of view, we are in particular interested in the diffusive substrate flux Ddc into the biofilm. dz z¼k Integrating (1) once gives
Jðcb ; kÞ :¼
dc lX 1 ðkÞ ¼ dz YD
Z
k
0
c dz; jþc
ð2Þ
where we choose our notation as to indicate that the flux depends on both, the bulk substrate concentration cb and the biofilm thickness k. We non-dimensionalize (1) using the scalings x ¼ kz ; s ¼ jc and sb ¼ cjb and obtain 2
d s 2
dx
¼ ak2
where a ¼
X 1l . YDj
s ; 1þs
sð1Þ ¼ sb ;
ds ð0Þ ¼ 0; dx
ð3Þ
Instead of (2) we obtain
jðsb ; kÞ :¼
Z 1 ds s 2 dx: ¼ ak dxx¼1 1þs 0
ð4Þ
Our aim is to find a closed form analytical approximation of (4).
Table 1 Model parameters from [16]. Symbol
l D X1
j Y
Parameter Maximum growth rate Diffusion coefficient Biomass density Monod constant Yield coefficient
Units 1
T L2 T1 M L3 M L3 –
Value 1.65 105 s1 1.0 105 cm 2 s1 2.5 103 g cm3 3.9 106 g cm3 0.071
1486
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
3. Analytical series solution 3.1. Derivation Starting point is (3), re-written in the form 2
2
d s
d s þ s 2 ak2 s ¼ 0; 2 dx dx
ds ð0Þ ¼ 0: dx
sð1Þ ¼ sb ;
ð5Þ
We define for p 2 [0, 1] the homotopy
"
v ðx; pÞ
H1 ðx; pÞ ¼ ð1 pÞ
# " # @2 @2 @2 @2 2 v ðx; pÞ v ðx; pÞ v ðx; 0Þ þ p v ðx; pÞ þ v ðx; pÞ v ðx; pÞ ak v ðx; pÞ ¼ 0; @x2 @x2 @x2 @x2
d v ð0; pÞ ¼ 0: dx
v ð1; pÞ ¼ sb ;
ð6Þ
It is easy to verify that for p = 1 this is equivalent with (3). Moreover, every function that satisfies the boundary conditions is a valid solution of (6) for p = 0. We expand v(x; p) in a Taylor series about p = 0 to obtain
v ðx; pÞ ¼ v ðx; 0Þ þ
þ1 X
v ðmÞ ðx; pÞjp¼0 m!
m¼0
!
pm :
ð7Þ
Formally, provided there is sufficient regularity in the problem and the series is convergent, we can write, therefore,
v ðx; 1Þ ¼ v ðx; 0Þ þ
þ1 X
v ðmÞ ðx; pÞjp¼0
m¼0
m!
ð8Þ
:
We call the terms under the sum the mth order deformation derivative, defined as
sm ðxÞ :¼
v ðmÞ ðx; pÞjp¼0 m!
ð9Þ
:
The solution of (5) is then
sðxÞ ¼ s0 þ
þ1 X
sm ðxÞ;
ð10Þ
m¼1
where s0(x) = v(x; 0). The material flux into the biofilm can be computed from this according to the first equality in (4) by differentiation and multiplication with a constant. The deformation derivatives can be computed in analogy to the theory summarized in [10]. For our specific problem we have the following algorithm. Theorem 3.1. For s0(x) = sb, the deformation derivatives sm associated with homotopy H1 are obtained recursively as solutions of the boundary value problems m X k¼0
2
sk
d smk 2
dx
2
þ
d sm1 2
dx
ak2 sm1 ¼ 0;
sm ð1Þ ¼ 0;
dsm ð0Þ ¼ 0: dx
ð11Þ
Proof. Substituting (7) into (6) and applying Leibniz’ Rule, we obtain for the left hand side of (6)
m X m L:H:S ¼ ð1 pÞðlÞ l l¼0
"
@2 @2 v ðx; pÞ 2 v ðx; pÞ v ðx; pÞ 2 v ðx; 0Þ @x @x
#ðmlÞ
or
L:H:S ¼
m X m k¼0
k
v ðkÞ ðx; pÞ
m1 Xm 1 @ 2 ðmkÞ @ 2 ðmk1Þ ðkÞ v ðx; pÞj m v ðx; pÞ v ðx; pÞjp¼0 : p¼0 @x2 @x2 k k¼0
For the right hand side we obtain
" #ðmlÞ m 2 X m @2 ðlÞ @ 2 R:H:S ¼ v ðx; pÞ þ v ðx; pÞ 2 v ðx; pÞ ak v ðx; pÞ ðpÞ @x2 @x l l¼0
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
1487
or
m1 X m1 @ 2 ðm1Þ @2 v ðx; pÞj m v ðkÞ ðx; pÞ 2 v ðm1kÞ ðx; pÞjp¼0 þ mak2 v ðm1Þ ðx; pÞjp¼0 : p¼0 2 @x @x k k¼0
R:H:S ¼ m
Rewriting this in terms of the deformation derivatives sm(x) this becomes
R:H:S ¼ m!
" 2 d sm1 2
dx
þ
m1 X
#
2
sk
d sm1k
ak2 sm1
2
dx
k¼0
and
" L:H:S ¼ m!
m X
2
d smk
sk
2
dx
k¼0
m1 X
2
d smk1
sk
2
dx
k¼0
# :
Thus, for the mth deformation derivative we obtain with the boundary conditions in Eq. (6) the linear boundary value problem m X
2
sk
k¼0
2
d smk dx
þ
2
d sm1 2
dx
ak2 sm1 ¼ 0;
dsm ð0Þ ¼ 0: dx
sm ð1Þ ¼ 0;
Note that for given s1 ; . . . ; sm1 , the function sm can be computed by integrating twice. Since we choose s0(x) = sb = const, the mth order derivative sm is an even polynomial of degree less than or equal 2m. Thus, s(x) is a conventional power series for the solutions of (3). Its coefficients can be computed conveniently with computer algebra software. We use Maple for this task. We now consider the homotopy
"
# " # @2 @2 @2 @2 2 H2 ðx; pÞ ¼ ð1 pÞ v ðx; pÞ 2 v ðx; 0Þ þ p 2 v ðx; pÞ þ v ðx; pÞ 2 v ðx; pÞ ak v ðx; pÞ ¼ 0; @x2 @x @x @x d v ð0; pÞ ¼ 0: dx
v ð1; pÞ ¼ 1;
ð12Þ
Theorem 3.2. For s0(x) = sb, the deformation derivatives sm associated with homotopy H1 are obtained recursively as solutions of the boundary value problems 2
d sm 2
dx
¼
m1 X
2
sk
k¼0
d
s 2 m1k
dx
þ ak2 sm1 ;
sm ð1Þ ¼ 0;
dsm ð0Þ ¼ 0: dx
Proof. As in the previous proof, we substitute (7) into (12) and use Leibniz’s rule to obtain
" #ðmlÞ m 2 X m @2 ðlÞ @ L:H:S ¼ v ðx; pÞ 2 v ðx; 0Þ ð1 pÞ @x2 @x l l¼0 or
L:H:S ¼
@ 2 ðmÞ @ 2 ðm1Þ v ðx; pÞj m v ðx; pÞjp¼0 : p¼0 @x2 @x2
Moreover, we obtain
" #ðmlÞ m 2 X m @ 2 v ðx; pÞ ðlÞ @ v ðx; pÞ 2 R:H:S ¼ þ v ðx; pÞ ak v ðx; pÞ ðpÞ @x2 @x2 l l¼0 or
R:H:S ¼ m
m1 X m1 @ 2 ðm1Þ @ 2 ðm1kÞ ðkÞ v ðx; pÞj m v ðx; pÞ v ðx; pÞjp¼0 þ mak2 v ðm1Þ ðx; pÞjp¼0 : p¼0 @x2 @x2 k k¼0
Re-writing this in terms of the deformation derivative, we obtain
" R:H:S ¼ m!
2
d
s 2 m1
dx
þ
m1 X k¼0
#
2
sk
d
s 2 m1k
dx
2
ak sm1
ð13Þ
1488
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
and
" L:H:S ¼ m!
# 2 d : s s 2 m 2 m1 dx dx 2
d
Thus, 2
d sm 2
dx
¼
m1 X k¼0
sk
d
2
dx
s 2 m1k
þ ak2 sm1 ;
sm ð1Þ ¼ 0;
dsm ð0Þ ¼ 0: dx
Again one obtains sm(x) as an even polynomial of degree less than or equal 2m. 3.2. Numerical comparison The recursions in Theorems 3.1 and 3.2 allow us to compute approximations to the solutions of (3) but do not give any insight how fast, if at all, these approximations converge. Every finite order approximation is a polynomial approximation of the solution. One problem that may arise is that the power series representation of the Monod nonlinearity in (3) has only a finite radius of convergence which indicates a finite radius of convergence for the polynomial approximation of the solution as well. This was accounted for to some extent by providing two homotopies, one to be used for sb < 1 and one for sb > 1. By construction, the coefficients of sm(x) according to Theorems 3.1 and 3.2 depend on the parameters ak2 and sb. Therefore, it P can be expected that these affect the convergence of the series sðxÞ ¼ 1 m¼0 sm ðxÞ. Due to the algebraic complexity, a rigorous convergence proof is difficult. Moreover, mere convergence of the series does not automatically imply that truncating the series after few terms leads to a good quantitative approximation of the solution of the underlying two-point boundary value problem. Therefore, we carry out a computational comparison of finite-term approximations according to (3.1) and (3.2) with a numerical solution. The latter is computed with Maple’s built in Finite Difference method. In fact, since our original goal was to determine the flux (4) of substrate into the reactive layer, we do not compare s(x) but the approximate flux s0 (1) against the numerically computed one. The set of model parameters in this study is taken from [16] and listed in Table 1. The lth order fluxes according to the analytical approximation are computed for k and cb = sbj as
J l ðcb ; kÞ ¼
l jX
k
s0m ð1Þ;
ð14Þ
m¼0
where s0m ð1Þ is calculated with the computer algebra software Maple according to the recursion in Theorem 3.1 for cb < j and according to Theorem 3.2 for cb > j. The results are plotted in Fig. 1 for four selected values of cb and for selected flux approximations Jl. In all cases we observe an acceptable agreement of Jl with the numerical solution J num ¼ dc ðkÞ for small dz k<~ k 2 105 cm only. For larger values of k, this approximation breaks down, with the possible exception of the largest bulk concentration value tested. From an application point of view the small range of k for which the homotopy based methods work renders the approximation unusable. Bacterial cells have a diameter of order of k = 104 cm, i.e. biofilms, as bacterial layers, are considerably thicker than ~ k. Although this is a negative result, we can use the non- or slow-convergent homotopy based approximations to motivate an easy to evaluate flux approximation in the next section. 4. Quadratic approximation to the series solutions Upon inspection of the deformation derivatives sm(x) for m 6 13 according to Theorem 3.1 and for i 6 10 according to Theorem 3.2, we propose a new approximation. It is based on the observation that sm(x) can be expressed as a quadratic polynomial in x plus Oðk4 Þ terms. We treat the cases of Theorem 3.1 (used for sb > 1) and Theorem 3.2 (used for sb < 1) separately but will see that both will lead to the same compact formula for a flux approximation. This will be compared against the numerical solution. 4.1. Derivation Lemma 4.1. The solution sm(x) of the generalized mth boundary value problem (11) for m P 1 has the form
sm ðxÞ ¼
ð1Þm1 2 2 ak x 1 þ Oðk4 Þ: 2sm1 b
Proof. We show the assertion for m = 1 explicitly and use mathematical induction for the general case. The boundary value problem (11) for m = 1 is
ð15Þ
1489
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494 8e-09
2e-10
6e-09
J [g/cm44]
1.5e-10
J [g/cm4]
numerical J2 J=9 J13
numerical J2 J=8 J10
1e-10
4e-09
2e-09
5e-11
0
0 0
1e-05
2e-05
3e-05
4e-05
0
5e-05
1e-05
lambda [cm]
6e-09
8e-10
2e-05
3e-05
4e-05
3e-05
4e-05
lambda [cm]
numerical J2 J=9 J13
numerical J2 J=8 J10
6e-10
J [g/cm4]
J [g/cm4]
4e-09
4e-10
2e-09 2e-10
0
0 0
1e-05
2e-05
3e-05
4e-05
5e-05
0
lambda [cm]
1e-05
2e-05
lambda [cm]
Fig. 1. Comparison of flux approximation according to (14) and Theorem 3.2 for sb = cb/j < 1 (left column) and Theorem 3.1 sb = cb/j > 1 (right column). Note: j = 3.9 106 g/cm3, see Table 1.
s0 s001 þ s1 s000 þ s000 ak2 s0 ¼ 0;
s1 ð1Þ ¼ 0;
s01 ð0Þ ¼ 0:
Since s0(x) = sb, simple integration gives
s1 ðxÞ ¼
ak2 2 ðx 1Þ; 2
which is of form (15). Assume that the assertion holds for l = 1, . . . , m 1. From (15) we obtain
s00m1 ðxÞ ¼
ð1Þm2 2 ak þ Oðk4 Þ: sm2 b
Substituting sm1 ðxÞ and s00m1 ðxÞ into (11) leads to the integration problem
s0 s00m ðxÞ þ s00m1 ðxÞ þ Oðk4 Þ ¼ 0;
sm ð1Þ ¼ 0;
s0m ð0Þ ¼ 0:
Thus
sm ðxÞ ¼
ð1Þm1 2 2 ak ðx 1Þ þ Oðk4 Þ: 2sm1 b
Proposition 4.2. The analytical approximate solution s(x) of the Monod boundary value problem (3) according to homotopy (6) for sb > 1 has the form
1490
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
sðxÞ ¼ sb þ
ak2 ðx2 1Þ sb þ Oðk4 Þ: 2 1 þ sb
ð16Þ
Moreover,
s0 ð1Þ ¼
sb ak2 þ Oðk4 Þ: 1 þ sb
ð17Þ
Proof. With the previous lemma we obtain
sðxÞ ¼ sb þ
1 ak2 ðx2 1Þ X ð1Þm1 þ Oðk4 Þ: 2 sm1 b m¼0
For sb > 1, this is
sðxÞ ¼ sb þ
ak2 x2 1 sb þ Oðk4 Þ: 2 1 þ sb
which is in the required form. For the derivative at k = 1 it follows
s0 ð1Þ ¼
sb ak2 þ Oðk4 Þ: 1 þ sb
Thus, by dropping the Oðk4 Þ terms in Eq. (16) a candidate for a closed form analytical solution is obtained for the bulk concentration sb > 1. We repeat the same analysis for the solution s(x) obtained with homotopy (12). Lemma 4.3. The solution sm(x) of the generalized mth boundary value problem (13) for m P 1 has the form
sm ðxÞ ¼
ak2 ð1 x2 Þ 4 ð1Þm sm b þ Oðk Þ: 2
ð18Þ
Proof Again we use mathematical induction. For m = 1 we obtain
s0 s000 þ s001 ak2 s0 þ Oðk4 Þ ¼ 0;
s1 ð1Þ ¼ sb ;
s01 ð0Þ ¼ 0
and thus
s1 ¼
ak2 ðx2 1Þ sb þ Oðk4 Þ: 2
This is the same expression obtained in the previous case. The assertion is assumed to hold for l = 1, 2, . . . , m 1, which means
sm1 ðxÞ ¼
ak2 ð1 x2 Þ þ Oðk4 Þ ð1Þm1 sm1 b 2
is true. The second order derivative is then
s00m1 ðxÞ ¼
ak2 þ Oðk4 Þ: ð1Þm sm1 b 2
Substituting sm1 and s00m1 into (13) we obtain
s00m ¼ s0 s00m1 ;
sm ð1Þ ¼ 0;
sm ð0Þ ¼ 0
from which
sm ðxÞ ¼
ak2 ð1 x2 Þ 4 ð1Þm sm b þ Oðk Þ: 2
Proposition 4.4. The analytical approximate solution s(x) of the Monod boundary value problem (3) according to homotopy (12) for sb < 1 has the form
sðxÞ ¼ sb þ
ak2 ðx2 1Þ sb þ Oðk4 Þ: 2 1 þ sb
ð19Þ
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
1491
Moreover,
s0 ð1Þ ¼
sb ak2 þ Oðk4 Þ: 1 þ sb
ð20Þ
Proof. Using sm(x) from (18) gives
sðxÞ ¼ sb þ
1 ak2 ð1 x2 Þ X 4 ð1Þm sm b þ Oðk Þ: 2 m¼1
Thus, for sb < 1,
sðxÞ ¼ sb þ
ak2 ðx2 1Þ sb þ Oðk4 Þ: 2 1 þ sb
The derivative at k = 1 is then obtained as
s0 ð1Þ ¼
sb ak2 þ Oðk4 Þ: 1 þ sb
Again, dropping the Oðk4 Þ terms in (19) we derive a candidate for a closed form approximation of the solution of the Monod boundary value problem. We point out that the approximation candidates obtained for Theorem 3.1 (for sb > 1) and Theorem 3.2 (for sb < 1) take the same expression, which is continuous at sb = 1. 4.2. Numerical verification We repeat the numerical verification for the flux approximation obtained in Propositions 4.2 and 4.4. After dropping Oðk4 Þ term and converting back into physical units, we obtain for the derivative of the substrate concentration at the boundary
J analytical ¼
lX 1 YD
k
cb
j þ cb
ð21Þ
:
We compare this approximation against the numerical solution of (1), which we calculated with a standard finite difference method on a very fine grid. Converting the analytical expression back into physical units, we obtain
J num ¼
dc ðkÞ; dz
where the derivative is approximated by a finite difference from the concentration values in the last two grid points. We use the parameters from [16], summarized in Table 1. The biofilm thickness is varied over the range 20 lm < k < 200 lm. The bulk substrate concentration ranges over 0 6 cb 6 2 104 g/cm3. These are realistic ranges for many biofilm applications. For three selected biofilm thicknesses, k = 74 lm, 110 lm, 164 lm, the results are plotted in Fig. 2. The analytical approximation is found to be in good agreement with the numerical solution, in particular for smaller k and large cb. The quality of the approximation decreases as k becomes large for small cb. Qualitatively the same results were obtained with a second set of parameters, taken from [19] (data not shown). In all cases the analytical approximation is slightly larger than the numerical solution. This is in accordance with the following theoretical result. Proposition 4.5. The analytical approximation of the flux according to (17) or (20) after dropping the Oðk4 Þ terms overestimates the flux of the exact solution of (3). Proof. We first show that the Oðk4 Þ approximation ~s according to (16),
~sðxÞ ¼ sb þ
ak2 ðx2 1Þ sb 2 1 þ sb
is a lower estimate for s(x), i.e. that ~sðxÞ 6 sðxÞ. Indeed, we have on the one hand
~s00 ¼ ak2
sb 1 þ sb
and on the other hand
ak2
1 1 þ 12 ak2 ðx2 1Þ 1þs ~sðxÞ sb b : ¼ ak2 s 2 1 1 þ ~sðxÞ 1 þ sb 1 þ 2 ak ðx2 1Þ b 2 ð1þsb Þ
1492
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
0.001
0.0008
J [g/cm4]
0.0006
0.0004
0.0002 lambda=0.0164cm lambda=0.0110cm lambda=0.0074cm
0 0
5e-05
0.0001
0.00015
0.0002
Cb [g/cm3] Fig. 2. Analytical vs numerical computation of the specific substrate flux J for three values of biofilm thickness k as a function of the bulk concentration cb: The analytical solution is drawn as solid line, the numerical solution is represented by marks.
Since
sb 1þsb
< 1, the numerator in the last expression is larger than the denominator, thus
~s00 6 ak2
~sðxÞ : 1 þ ~sðxÞ
Hence, by standard comparison theorems for two-point boundary value problems, e.g. as in [2], we have
~sðxÞ 6 sðxÞ: With ~sð1Þ ¼ sð1Þ ¼ sb , it follows that
~s0 ð1Þ P s0 ð1Þ:
5. Simple applications In order to demonstrate how an analytical flux approximation can simplify biofilm models we include two examples. 5.1. Stability of biofilm growth The standard single-species/single substrate biofilm model [19] can be written in the form
dk ¼ dt
Z 0
k
lX 1
c kl X 1 dz dðkÞk; jþc
ð22Þ
where kl is the cell death rate and d(k) the rate of biofilm detachment. The substrate concentration c is given by (1). Together (1) and (22) form a free boundary value problem for c. From an application point of view, however, usually one is more interested either in the biofilm thickness k or in the amount of substrate converted by the biofilm, i.e. the flux of substrate into the biofilm. Using (21) and sb ¼ cjb , the problem is simplified to the ordinary differential equation
dk cb k ðkl þ dðkÞk ¼ FðkÞ; ¼A dt j þ cb
ð23Þ
where the constant A contains the original model parameters and does not depend on cb and k. It represents the growth rate of the biofilm and accounts for diffusion of the substrate in the biofilm. 0 b We compute F 0 ðkÞ ¼ jAc kl dðkÞ d ðkÞk. The stability of the model, therefore, depends on the parameters and the speþcb cific detachment rate function. The most commonly used detachment function is the so-called L2 detachment model, in which case d(k) = dk [19]. In this case the biofilm model (22), (1) becomes with our flux approximation the logistic equation
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
dk ¼ kðB dkÞ; dt
B :¼
Acb
j þ cb
1493
kl :
If the compounded net growth rate B is positive, a biofilm will establish, i.e. the equilibrium k = 0 is unstable and the equilibrium k = B/d is asymptotically stable. We note that such a logistic growth model was used ad hoc to be fitted against experimental data for detachment experiments in [17]. Our results give an interpretation of the coefficient B that allow us to relate it to the bulk substrate concentration and hence to environmental conditions. It is well understood and was demonstrated in experiments [13,17] that in fact the detachment rate depends on the reactor hydrodynamics. More specifically, it can be assumed that the detachment rate d(k) depends on the shear rate, and thus on the reactor type. A frequently used laboratory reactor is the so-called roto-torque reactor, where biofilm grows between two concentric cylinders, one of which rotates. This can be approximated by biofilm growing between two parallel plates, a distance D apart. One of the plates is stationary, the other one is assumed to move at a constant velocity u. As long as u is small enough that no Couette–Taylor instabilities occur, the flow profile is linear with slope u/(D k). Thus, the detachment rate is in this case proportional to u/(D k) [20]. Hence, from (22), (1) for the roto torque reactor with shear rate dependent detachment rate we have
dk c ¼k B ¼: FðkÞ; dt Dk where c is proportional to the reactor’s rotation speed. We have
F 0 ðkÞ ¼ B
ck c : D k ðD kÞ2
The trivial equilibrium k = 0 (no biofilm, washout) is stable, i.e. F0 (k) < 0 if the flow rate, i.e. c, is fast enough, or if the growth rate, i.e. B, is small enough. This is the case if cb is small. Depending on parameters, a second, positive, equilibrium k⁄ exists if 0 < k :¼ D Bc, i.e. if the trivial equilibrium is unstable. Here we can assume that B > 0 because otherwise no biofilm growth would be possible, even in the absence of detachment. If it exists, k⁄ is asymptotically stable. 5.2. Reactor scale analysis Similar as [13], we consider a vessel of volume V that is filled with biofilm carrier chips that provide additional surface area for biofilms to form. This mimics a Moving Bed Biofilm Reactor of the type used in wastewater engineering. The inflow and outflow rate to and from the vessel is Q. On inflow, the substrate concentration is c1. In the vessel, we assume that the substrate is completely mixed. The system is modeled by (22) and (1), together with a reactor mass balance for the substrate concentration in the bulk cb
dcb Q dc ¼ ðc1 cb Þ E ðkÞ; V dz dt where the new parameter E subsumes reaction–diffusion and reactor parameters. The last term models the consumption of substrate by the biofilm. Using (21) as before, the model reduces to the two-dimensional ordinary differential equation
dk cb k ¼A ðkl þ dðkÞÞk; dt j þ cb dcb cb k ; ¼ aðc1 cb Þ b dt j þ cb
ð24Þ ð25Þ
where parameters a, b subsume the other model parameters. This model can be studied in the phase plane with elementary techniques after the detachment rate function d(k) is specified. For example, for the standard choice d(k) = dk as above, the Jacobian of this system is obtained as
2 F ðk; cb Þ ¼ 4 0
Acb
jþcb
kl 2dk jbcþcb
b
3
Ajk ðjþcb Þ2
j a ðjbk þc Þ2
5:
b
c1 The system permits a trivial equilibrium k ; cb ¼ ð0; c1 Þ, for which we obtain the eigenvalues k1 ¼ A jþc c and 1 k2 ¼ a < 0. Thus, this equilibrium is unstable, i.e. a biofilm will be able to establish itself, if the bulk substrate concentration c1 c1 is high enough that the maximum growth rate is bigger than the cell death rate, i.e. A jþc > c, otherwise it is stable. In the 1 ⁄⁄ ⁄⁄ former case, a positive non-trivial equilibrium (k , c ) is obtained as the unique intersection of the strictly increasing funccb tion ^ k ¼ Ad jþc cd with the strictly decreasing function ~ k ¼ ab ðC 1 cb Þðj þ cb Þ. Then k⁄⁄ > 0 and c b < c1 . The Jacobian in this b case becomes 0
F k
; c b
2 ¼4
dk
3
Ajk ðjþc Þ2 b
j
bk Ab ðc þ dk Þ a ðjþc Þ2 b
5
1494
F. Abbas, H.J. Eberl / Applied Mathematics and Computation 218 (2011) 1484–1494
with negative trace and positive determinant. Thus k ; c is asymptotically stable. b Hence, in all three applications that we briefly sketched here, the analytical flux approximation allows us to reduce a hybrid ODE/PDE model to a simple ordinary differential equation that can be studied with elementary techniques. 6. Conclusion Two different homotopy equations (6) and (12) were proposed to obtain series solutions for the Monod two point boundary value problem. Both lead to approximations that are irrelevant for practical purposes. However, using that the biofilm thickness k as a small parameter enters these approximations to the 4th power, we were able to obtain an easy to evaluate Oðk4 Þ analytical approximation of the infinite order approximations, which was shown to be an accurate approximation of the numerically computed substrate flux into the biofilm. This approximation will allow to simplify modeling tasks, such as upscaling of biofilm models and stability analysis, by reducing model complexity. Acknowledgement This research was supported in parts by the Higher Education Commission of Pakistan through a doctoral scholarship for FA. References [1] S. Abbasbandy, The application of the homotopy analysis method to nonlinear equations arising in heat transfer, Phys. Lett. A 360 (2006) 109–113. [2] S. Carl, S. Heikkilä, Nonlinear Differential Equations in Ordered Spaces, Chapman & Hall, 2000. [3] M.S.H. Chowdhury, I. Hashim, Analytical solutions to heat transfer equations by homotopy-perturbation method revisited, Phys. Lett. A 372 (2008) 1240–1243. [4] M. Duman, Asymptotic expansions for the Sturm–Liouville problem by homotopy perturbation method, Appl. Math. Comput. 216 (2010) 492–496. [5] D.D. Ganji, The application of He’s homotopy perturbation method to nonlinear equations arising in heat transfer, Phys. Lett. A 355 (2006) 337–341. [6] J.H. He, Comparison of homotopy perturbation method and homotopy analysis method, Appl. Math. Comput. 156 (6) (2004) 527–539. [7] J.H. He, Homotopy perturbation method for solving boundary value problems, Phys. Lett. A 350 (2006) 87–88. [8] S. Liang, D.J. Jeffrey, Comparison of homotopy analysis method and homotopy perturbation method through an evolution equation, Commun. Nonlinear Sci. Numer. Simul. 14 (12) (2009) 4057–4064. [9] S.J. Liao, Beyond Perturbation: Introduction to The Homotopy Analysis Method, Chapman and Hall, 2004. [10] S.J. Liao, Notes on the homotopy analysis method: some definitions and theorems, Commun. Nonlinear Sci. Numer. Simul. 14 (4) (2009) 983–997. [11] S.J. Liao, An optimal homotopy-analysis approach for strongly nonlinear differential equations, Commun. Nonlinear Sci. Numer. Simul. 15 (8) (2010) 2003–2016. [12] S.J. Liao, On the relationship between the homotopy analysis method and Euler transfer, Commun. Nonlinear Sci. Numer. Simul. 15 (6) (2010) 1421– 1431. [13] A. Masic, J. Bengtsson, M. Christensson, Measuring and modeling the oxygen profile in a nitrifying Moving Bed Biofilm Reactor, Math. Biosci. 227 (1) (2010) 1–11. [14] A. Rafiq, M. Ahmed, S.A. Hussain, A general approach to specific second order ordinary differential equations using homotopy perturbation method, Phys. Lett. A 372 (2008) 4973–4976. [15] B.E. Rittmann, P.L. McCarty, Model of steady state biofilm kinetics, Biotechnol. Bioeng. 22 (11) (1980) 2343–2357. [16] C.T. Skowlund, Effects of biofilm growth on a steady state biofilm model, Biotechnol. Bioeng. 35 (5) (1990) 502–510. [17] Y.P. Tsai, Impact of flow velocity on the dynamic behaviour of biofilm bacteria, Biofouling 21 (5) (2005) 267–277. [18] Y.G. Wang, H.F. Song, D. Li, Solving two-point boundary value problems using combined homotopy perturbation method and Green’s function method, Appl. Math. Comput. 212 (2009) 366–376. [19] O. Wanner, H. Eberl, E. Morgenroth, D. Noguera, C. Picioreanu, B. Rittmann, M.C.M. van Loosdrecht, Mathematical Modeling of Biofilms, IWA Publishing, 2006. [20] T. Xu, Numerical and statistical analysis of detachment forces in biofilm colonies, M.Sc. Thesis, Dept. Mathematics and Statistics, University of Guelph, 2005.