Net growth rate of continuum heterogeneous biofilms with inhibition kinetics.

Biofilm systems can be modeled using a variety of analytical and numerical approaches, usually by making simplifying assumptions regarding biofilm heterogeneity and activity as well as effective diffusivity. Inhibition kinetics, albeit common in experimental systems, are rarely considered and analytical approaches are either lacking or consider effective diffusivity of the substrate and the biofilm density to remain constant. To address this obvious knowledge gap an analytical procedure to estimate the effectiveness factor (dimensionless substrate mass flux at the biofilm-fluid interface) was developed for a continuum heterogeneous biofilm with multiple limiting-substrate Monod kinetics to different types of inhibition kinetics. The simple perturbation technique, previously validated to quantify biofilm activity, was applied to systems where either the substrate or the inhibitor is the limiting component, and cases where the inhibitor is a reaction product or the substrate also acts as the inhibitor. Explicit analytical equations are presented for the effectiveness factor estimation and, therefore, the calculation of biomass growth rate or limiting substrate/inhibitor consumption rate, for a given biofilm thickness. The robustness of the new biofilm model was tested using kinetic parameters experimentally determined for the growth of Pseudomonas putida CCRC 14365 on phenol. Several additional cases have been analyzed, including examples where the effectiveness factor can reach values greater than unity, characteristic of systems with inhibition kinetics. Criteria to establish when the effectiveness factor can reach values greater than unity in each of the cases studied are also presented.

Introduction

The successful design of large-scale bioreactors for pollutant degradation requires the ability to accurately predict the growth of microorganisms. A great number of theoretical investigations of the growth rate of one-species or mixed species biofilms have focused on the difficulty of transferring substrates and products between the fluid phase and the cells inside a biofilm. These analyses were based on the theory of mass transport and diffusion in porous media but rarely considered inhibition. In contrast, many experimental investigations have reported inhibition kinetics of different types, e.g., substrate inhibition,[1-4] product inhibition,[5-7] substrate and product inhibition,[8,9] inhibition by other compounds present in the fluid phase,[10] inhibition effects of heavy metals ions,[11] and wastewater treatment under salt-affected conditions,[12] among others. Equations were usually fitted to experimental data in all of these field or laboratory studies, without developing predictive scenarios of what would happen in different situations (for example, as a result of product inhibition), and only some specifically addressed biofilm systems.[5,11] Modeling studies on substrate utilization and inhibition effects in biofilms are few and far between.[13-16] It is, therefore, of interest to model the effect of the coupled processes of substrate and product diffusion on the net rate of biofilm growth when different types of inhibition mechanisms are considered.

The objectives in the present study were to extend a procedure to estimate analytically the effectiveness factor (dimensionless substrate mass flux at the biofilm-fluid interface) for a continuum heterogeneous biofilm with multiple limiting substrates Monod kinetics[17] to different types of inhibition kinetics. We tested a variety of scenarios when (i) the substrate is the limiting component, (ii) the inhibitor is the limiting species, or (iii) the substrate is also the inhibitor. The model describes a heterogeneous biofilm with variable distribution of biofilm density, activity, and effective diffusivity.

Model development

Inhibition modeling

The effect of an inhibitory component on the specific growth rate, r, is given by:[18,19]where qmax is the maximum specific growth rate, CA is the concentration of the key substrate A, KA is the Monod half rate constant for substrate A, CI is the concentration of the inhibitor I, and KI is the concentration of the inhibitor resulting in 50% inhibition of the maximum rate (all variables are described in detail in Box 1).

Taking into account that both the effective diffusivities of substrate and inhibitory component and the biofilm density vary with the position (x) in the biofilm, the mass balances for the substrates and inhibitors, at steady state (pseudo-steady state), are given by:where DfA and DfI are the effective diffusivities of the substrate and inhibitor, Xf is the biofilm density, and YA and YI are the biomass yield coefficients of the key substrate A and inhibitor I, respectively.

The following assumptions are made in Eqs. (2) and (3):

The biofilm is a continuum.

Substrate and inhibitor are transferred by diffusion only, in accordance with Fick’s law.[20]

Microorganisms consume the substrate and inhibitor at a rate according to Double-Monod inhibition kinetics.

Biofilm density and substrate effective diffusivity change in the x direction.

Steady state (pseudo-steady state) conditions apply and the rate of substrate consumption is higher compared to the rate of biofilm growth.

External resistances to mass transfer are neglected. Therefore, the substrate and inhibitor concentrations at the biofilm-fluid interphase, CAs and CIs, respectively, are known.

The substratum where the biofilm grows is impermeable.

Under these assumptions, the appropriate boundary conditions for Eqs. (2) and (3) are:with Lf representing the average biofilm thickness.

The concentrations of the substrate, A, and inhibitory component, I, vary as a function of biofilm depth in the continuous heterogeneous biofilm (Fig. 1). The DfA(x), DfI(x), and Xf(x) profiles and dimensionless differential equations describing the system can be found in Supplementary Information 1.

Fig. 1

Schematic representation of the continuum heterogeneous biofilm model with concentration profiles for substrate A and inhibitory component I. CAs concentration of substrate (A) at the surface of the biofilm, kg/m3, CIs concentration of inhibitor (I) at the surface of the biofilm, kg/m3, Lf biofilm thickness. In black: impermeable substratum

Case study

We studied three possible scenarios:

Case (a): If the substrate (A) is the limiting component rather than the inhibitor, which can be a reaction product, thenwithTherefore, Eq. (S1–14) results in (see Supplementary Information 1):and

Case (b): If the inhibitor is the limiting substrate, rather than the substrate, thenwithTherefore, Eq. (S1–14) results inand

Case (c): This is a special situation of case (b), where there is a single substrate that is also the inhibitor.The dimensionless kinetic rate expression in this case is (see Eq. (S1–14)):with

Mass balance differential equations

The mass balance differential equations for the corresponding limiting component (in each case), considering valid Eq. (S1–13), are:

Case (a):

Case (b):

Case (c):

Particular conditions

Analyzing the asymptotic solutions of the different cases studied for ϕ ≪ 1 (see Supplementary Information 2), there is a condition under which the first derivative of r*, evaluated at  = 1, is negative. In this situation the value of parameter σ (σA, σi), will be negative. This is a necessary but not sufficient condition to assure that there will be a region of the η versus ϕ curve, in which η values will be greater than one.

These particular situations are met in the following scenarios below.

Case (a): Substrate A is the limiting substrate and the inhibitor is another reactant, as in Eq. (9). The effectiveness factor, η (the ratio between the diffusion-limited substrate consumption rate and the substrate consumption rate that is not limited by diffusion), can reach values greater than unity when:

Case (b): The inhibitor is the limiting component, as in Eq. (13). Then:

Case (c): The same component is both substrate and inhibitor, as in Eq. (17).

In all these cases, the matching expressions (S2–18), (S2–22) and (S2–26) are able to predict the true value of the effectiveness factor. When η is greater than one, the diffusion limited substrate consumption rate is greater than the rate without diffusion resistance.

Results and discussion

Model setup

The continuum heterogeneous model was applied to different inhibition scenarios:

Case (a): The substrate A is the limiting substrate and the inhibitor is another reactant, which can be considered in two different ways:

Case (a1): Double Monod inhibition kinetics where substrate A is the limiting component and substrate I the inhibitor.

Case (a2): Double Monod inhibition kinetics where substrate A is the limiting component and the inhibitor I is a product.

Case (b): Double Monod inhibition kinetics where the substrate I is the limiting component and reacts with substrate A.

Case (c): Monod inhibition kinetics where the single component A is both substrate and the inhibitor.

Each case was solved under conditions (parameter values) where effectiveness factor values greater than one could be obtained, except for case (a2), which does not have this possibility (see below and Fig. 2b). The effectiveness factor as a function of the Thiele modulus (which is the ratio of the time scale of reaction to the time scale of diffusion) for case (a1) considering a continuum heterogeneous biofilm, η, a homogeneous biofilm, η0, and a continuum heterogeneous biofilm without inhibition, ηNI, was considered (Fig. 2a). The possibility of finding values of η greater than one is predicted by the condition established by Eq. (21):

Fig. 2

Effectiveness factor for double Monod kinetics with inhibition. Ψ = 0.5, κ = 4. a Case (a1): Substrate A is the limiting component; βA = 0.05, βI = 1.25, ΓI = 0.7. b Case (a2): The product acts as inhibitor βA = 0.1, βBI = 1, ΓB = −3. c Case (b): The inhibition compound is the limiting component; βA = 0.05, βI = 1.25, ΓA = 0.7. d Case (c): The same component is substrate and inhibitor; βA = 0.01, βAI = 1. η: effectiveness factor for a continuum heterogeneous biofilm with inhibition; ηo: effectiveness factor for a homogeneous system with inhibition; ηNI: effectiveness factor for heterogeneous biofilm without inhibition (single Monod kinetics)

For case (a2) the inhibitor is a process product. In this case, η values greater than one cannot be found because will always be positive, since is always negative, due to the fact that νI will be negative (Fig. 2b).

For Double Monod inhibition kinetics where the inhibitor is the limiting component, as in case (b), it is observed that in a continuum heterogeneous biofilm with Thiele modulus values up to 1, the effectiveness factor has values greater than one (maximum ≈ 1.1) (Fig. 2c). This case also obeys the criterion that establishes the condition to find η > 1.Therefore, 0.8 > 0.033.

Finally, Fig. 2d shows the effectiveness factor for case (c), Double Monod inhibition kinetics with a single substrate that is also the inhibitor. The parameters values were chosen in such a way to meet conditions that allow the possibility to find values of η > 1. Values of η as high as 1.1 were found. For comparison purposes only, the effectiveness factors for a homogeneous biofilm, η0, and a continuum heterogeneous biofilm without inhibition, ηNI, (Single Monod kinetics) are shown. As expected, the criterion to find η values greater than one is also obeyed.

Experimental validation

To assess the robustness of the new biofilm model with inhibition Monod kinetics, a test case was chosen where the substrate acts as nutrient as well as inhibitor (Table 1). The kinetic parameters were taken from the work of Chung et al.[21] who studied the growth of Pseudomonas putida CCRC 14365 during the biodegradation of phenol. The authors used the Haldane inhibition kinetics expression, which is a particular situation of the more general Monod kinetics. It therefore coincides with the Monod kinetics used in the present study, since the ratio between the substrate saturation constant and the inhibition constant (KP/KPI) is small (≈0.08).[11] The Monod growth rate with inhibition is given by Eq. (1). Developing the equation further, it is found thatWhen (KP/KPI) ≪ 1, then Eq. (24) becomes

Table 1

Parameters used for the test case of a continuum heterogenous biofilm model of Pseudomonas putida growth on phenola

Model Parameter	Designation	Value
Maximum specific growth rate	qmax	0.38 h−1
Phenolb saturation constant	K P	18.3 g of P m−3
Phenol inhibition constant	K PI	214.5 g of P m−3
Biomass yield	Y P	0.627 g of biomass g P−1
Physical properties
Phenol average diffusion coefficient	D¯_fP	2.7 × 10−6 m2 h−1
Average value of biomass density in biofilm	X¯_f	10,000 g of biomass/m3 of biofilm
Biofilm thickness	L f	75, 150 or 300 μm
Biofilm heterogeneity	Ψ	0.5
Relation DwP/αP	κ	4

Kinetic model:

Microbial kinetics at 30 °C

aExperimental data taken from Chung et al.[21]

bPhenol abbreviated to P

This Eq. (25) is the Haldane equation for substrate inhibition. In contrast to our study, Chung et al.[21] assumed a biofilm with a constant average density, , and a constant average substrate effective diffusivity, . In the present continuum heterogeneous biofilm model, the dimensionless effective diffusivity of phenol and biofilm density change with biofilm position () as:Therefore, the average biofilm density, , and the average substrate effective diffusivity, , are those indicated in Table 1.

Effect of biofilm thickness

Three cases were considered with varying biofilm thicknesses of 75, 150 and 300 μm. All the equations needed for solving the cases are given in Table 2. The net rate of biomass formation was predicted by the continuous heterogeneous biofilm model as a function of phenol concentration, CPs, at the biofilm-fluid interphase (Fig. 3a). The effectiveness factor, net rate of phenol consumption, and biomass production for phenol bulk concentrations in the range of 25 to 600 g/m3 were calculated using the values of the fundamental parameters given in Table 2. Initially, the biomass production rate increased very fast with increasing phenol concentration until a maximum was reached followed by a slow decrease (Fig. 3a). This behavior is in agreement with experimental findings.[1,21] The maximum rate was reached at lower phenol concentrations as the biofilm thickness decreased. However, for high phenol concentrations the thicker biofilm produced more biomass than the thinner one. The increase in the rate at a phenol concentration, CPs < 50 g/m3, was due to the fact that the rate tends towards a global first order reaction with respect to phenol concentration. The inhibition effect ([KPI /(KPI + CPs)] ≈ 1) was very small. However, when the phenol concentration was high (CPs greater than the substrate inhibition constant) the rate tended toward a negative reaction order (−1) with respect to phenol concentration. In other words, the reaction rate decreased as the substrate concentration increased. For phenol concentrations greater than 200 g/m3, inhibition played an important role in the overall rate of phenol consumption.

Table 2

Test case conditions

β_P=K_PC_Psβ_PI=K_PIC_Ps

r_PIs=q_maxX¯_fY_PC_PsK_P+C_PsK_PIK_PI+C_Ps

rPI*=β_P+1β_PI+1CP*β_P+CP*1β_PI+CP*

rPI*′1=β_Pβ_P+1-1β_PI+1ϕ^2=Lf2D¯_fPC_Psr_IPs

I_PI=β_P+1β_PI+1∫01CP*β_P+CP*1β_PI+CP*dCP*

ρ_PI=2DfP*1Xf*1I_PI^1∕2ϕ^*=ϕρ_PI

σ_PI=bψ^2X¯f2c0.2218rPI*′1FσPI*=σ_PIσPI2

d_PI=1-2σPI*

Data corresponding to:

aLf = 75 μm,

bLf = 150 μm,

cLf = 300 μm

Data without superscripts are the same for the three values of Lf

rIPs (g P/m3 h): rate of phenol consumption evaluated at the biofilm—fluid interface

rob(P) (g P/m3 min): net rate of phenol consumption in the biofilm

rob(B) (g B/m3 min): net rate of biomass production in the biofilm

Fig. 3

Test case: Model simulation when the substrate is also the inhibitor. Experimental data taken from Chung et al.[21] a Net rate of biomass production, rob, and b effectiveness factor, η, as a function of substrate biofilm-surface concentration, CPS, for different biofilm thicknesses, Lf, of 75, 150 and 300 μm

The fact that the maximum growth rate was reached at lower phenol concentrations as the biofilm thickness decreased is directly related to the internal diffusion effect. This is evident from the effectiveness factor observed at different phenol concentrations and for the three biofilm thicknesses (Fig. 3b). As can be seen, η is almost equal to unity for a 75 μm-thick biofilm. Here the effect of internal diffusion on the net rate of phenol consumption is at a minimum. For low substrate concentrations and a thin biofilm, both effects (internal diffusion resistance and inhibition) are negligible. However, as the biofilm thickness increases, internal diffusion resistance becomes more and more important. At high substrate concentrations (greater than 200 g/m3), the inhibition effect becomes more important than the diffusion resistance effect, even for thick biofilms. Under these conditions, substrate diffusion resistance and inhibition negatively affect the rate of reaction and the net rate decreases with increasing CPs.

Finally, it is interesting to note that the effectiveness factor takes on values greater than one, for the three cases where CPs > 200 g/m. To assure that η will take on values greater than one (case (c)) requires:For example, when applying this condition to CPs = 300 g/m3, βPI = 0.715 and βP = 0.061, for any biofilm thickness one obtains:

Another interesting observation regarding the behavior of η at different CPs (Fig. 3b) is the following: for each biofilm thickness, the point where η = 1 corresponds to a CPs value, where the inhibition effect is exactly compensated for by the substrate diffusion resistance.

Effect of mass transport resistance

The substrate mass transport inside a biofilm plays an important role in the estimation of the net rate microbial growth rate. A simple analytical method to calculate this rate using a model that considers variable diffusion and activity in the biofilm is highly useful.

Recently, Connolly et al.[16] demonstrated the influence of transport steps in biofilm growth by experimentally determining the net growth rate of Escherichia coli MJK2 to estimate the specific reaction rate as the net ureolysis rate per unit biofilm volume. In this case a homogeneous biofilm model was assumed and the convective external transport and diffusion rate into the biofilm were accounted for using dimensionless characteristic numbers like Damkôhler (Da) and Peclet (Pe) for external mass transport and the Thiele modulus ϕ (Eqs. S1–7) for diffusion into the biofilm, respectively. The reported average Da and Pe numbers of 4 and 104, respectively, and the Thiele modulus of less than one suggested that urea hydrolysis was not strongly limited by either external or internal diffusion. A Michaelis–Menten rate expression was used for fitting the experimental data with the finite element method. The values of the urea kinetic expression were varied until the difference between the modeled and experimental urea fluxes was minimized. However, simple first order kinetics with respect to urea concentration was found to fit the data best. Using our semi-analytical procedure to calculate the effectiveness factor one obtains the intrinsic rate expression in a more straightforward manner. Considering Eqs. (S2–18) we calculated the effectiveness factor, taking into account the Thiele modulus values of Connolly et al.[16] (Table 1), and the value of parameter ρ (Eqs. (S2–10)) for this case, which was used for normalizing the Thiele modulus. Further, the corresponding values of η according to our procedure for minima, maxima and average values of ϕ were 0.999, 0.947, and 0.998, respectively, indicating little to no internal biofilm diffusion limitation. Hence our procedure is able to estimate the effectiveness factor in experimental studies and validate conclusions on internal diffusion effects in biofilms.

Global applicability of continuum heterogeneous biofilm model to substrate inhibition

Olivieri et al.[15] considered a well-mixed three-phase bioreactor where the growth of a granular- Pseudomonas sp. OX1 biofilm was supported by the oxidation of phenol. The authors analyzed the system with either oxygen or phenol (the inhibitor) as limiting substrate. As in the present model, the empirical correlation between the substrate relative diffusivities and biomass concentration by Fan et al.[22] was used and both the steady-state and local dynamic behavior of the system were characterized. However, as in other earlier studies a homogeneous biofilm was assumed, that is, the effective diffusivity of the substrate and the biofilm density were kept constant. The model dimensionless differential equations were solved numerically under steady-state conditions. Although the conditions in that study and the present model system were quite different, effectiveness factor and net biofilm growth rate had a similar relationship with biofilm thickness.

Similarly, Pseudomonas stutzeri OX1 in an airlift biofilm reactor displayed specific growth rates as a function of the initial phenol concentration that were very close to those in our test case, with a maximum growth rate at about 200 g phenol/m3.[14]

When modeling batch fermentation kinetics for succinic acid production by Mannheimia succiniciproducens, both substrate and product inhibition were considered using glucose as substrate.[23] The growth of the microorganism was expressed by a simplified Monod model since the ratio of substrate saturation and inhibition constants was 0.0127 ≪ 1 and hence the Haldane model could be applied. The inhibition by the product was taken into account in the growth rate (1 – CPro/CProCrit), as suggested by Levenspiel,[24] where CPro is the product concentration and CProCrit is the critical product concentration at which cell growth stops. The effects of initial glucose concentration on the growth rate of M. succiniciproducens followed a profile that would also be obtained if our procedure were applied.

In conclusion, the continuum heterogeneous biofilm model successfully predicted situations where Double Monod inhibition kinetics affects the biofilm growth rate. Analytical approximate solutions can account for limiting concentrations of substrate or inhibitor and biofilm thickness. The procedure does not require numerical simulation and calculations can be performed in a basic software package. The approach has global utility for biofilm systems of different scales ranging from microfluidic flow cells to bioreactors used in laboratories and full-scale applications where inhibition kinetics is frequently encountered. The analytical procedure is accessible to researchers from different disciplines, especially those experimenting with flow cells where substrate limitation and inhibition kinetics may exert considerable effects on biological phenomena.

Methods

The model was designed as described in Model Development and the supplementary Information, and the experimental data used to validate the model were taken from ref. 21.

Data availability

The authors declare that data supporting the findings of this study are available within the paper and its Supplementary Information.

Nomenclature

dimensionless substrate (i) concentration. Eq. (S1–4)

Ci concentration of substrate (i), kg/m3

Cisconcentration of substrate (i) at the surface of the biofilm, kg/m3

dPI defined by Eq. (S2–28)

Dfi surface average effective diffusivity of substrate (i), m2/s

dimensionless relative effective diffusivity of substrate (i), defined by Eq. (S1–4)

average effective diffusivity of substrate (i) in the biofilm, m2/s

average effective diffusivity of inhibitor (I) in the biofilm, m2/s

Dwi diffusivity of substrate (i) in the liquid medium, m2/s

IPI defined by Eq. (S2–17)

Ki Monod half rate constant for substrate (i), kg/m3

KI Inhibitor concentration giving 50% inhibition rate, kg/m3

Lf average biofilm thickness, m

qmax maximum specific growth rate, s−1

r specific growth rate, s−1

rbvolumetric growth rate, g biomass/(m3biofilm ·h)

rs reference reaction rate defined by Eq. (S1–9)

r* dimensionless rate of reaction defined by Eq. (S1–4) or (S1–14)

dimensionless rate of reaction defined by Eq. (8)

dimensionless rate of reaction defined by Eq. (12)

dimensionless rate of reaction defined by Eq. (15)

first derivative of with respect to at

first derivative of with respect to at

first derivative of with respect to at

riob average rate of substrate (i) consumption of the whole biofilm, kg/s m3

Xf biofilm density, kg/m3

average biofilm density along the (x) direction, kg/m3

dimensionless relative density defined by Eq. (S1–4)

xdistance from the bottom of the biofilm, m

x*dimensionless distance defined by Eq. (S1–4)

Yi yield coefficient for substrate (i), (kg microorganism/kg nutrient)

Greek letters

αi effective diffusivity of substrate (i), at the bottom of the biofilm, m2/s

βi dimensionless parameter for substrate (i) defined by Eq. (S1–4)

ϕ Thiele modulus (the ratio between a reference reaction rate in a homogeneous biofilm, which is not diffusion limited and has a density equal to the average value in the biofilm, and the diffusion rate) defined by Eq. (S1–7)

ϕ*normalized Thiele modulus, defined by Eqs. (S2–19), (S2–23) or (S2–27)

ΓI parameter defined by Eq. (S1–8)

ΓA parameter defined by Eq. (10)

η effectiveness factor (the ratio between the diffusion-limited substrate consumption rate and the substrate consumption rate that is not limited by diffusion) for a continuum heterogeneous biofilm

η0 effectiveness factor for a homogeneous biofilm

ηNI effectiveness factor for continuum heterogeneous biofilm without inhibition

κ parameter defined by Eq. (S1–17)

νi ratio between the substrate (i) yield coefficient and that of the limiting substrate

ρAIparameter defined by Eq. (S2–16)

σA parameter defined by Eq. (S2–2)

σI parameter defined by Eq. (S2–5)

Ψparameter defined by Eq. (S1–16)

Sub indexes

i for substrate A or inhibitory substrate I

sbiofilm surface conditions

I for the inhibitory component

P for substrate phenol

Supplementary information