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][6][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][14][15][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 q max is the maximum specific growth rate, C A is the concentration of the key substrate A, K A is the Monod half rate constant for substrate A, C I is the concentration of the inhibitor I, and K I 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 D fA and D fI are the effective diffusivities of the substrate and inhibitor, X f is the biofilm density, and Y A and Y I are the biomass yield coefficients of the key substrate A and inhibitor I, respectively.
The following assumptions are made in Eqs. (2) and (3) with L f 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 D fA (x), D fI (x), and X f (x) profiles and dimensionless differential equations describing the system can be found in Supplementary Information 1.

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, then with Therefore, Eq. (S1-14) results in (see Supplementary Information 1): and Case (b): If the inhibitor is the limiting substrate, rather than the substrate, then Therefore, Eq. (S1-14) results in and 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)): Mass balance differential equations. The mass balance differential equations for the corresponding limiting component (in each case), considering valid Eq. (S1-13), are: Particular conditions Analyzing the asymptotic solutions of the different cases studied for ϕ 2 ≪ 1 (see Supplementary Information 2), there is a condition under which the first derivative of r*, evaluated at C Ã A = 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.  Each case was solved under conditions (parameter values) where effectiveness factor values greater than one could be obtained, except for case (a 2 ), 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 (a 1 ) 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):

RESULTS AND DISCUSSION
in this case 0:7 2:25 ¼ 0:311> 0:05 1:05 For case (a 2 ) the inhibitor is a process product. In this case, η values greater than one cannot be found because r Ã0 A ð1Þ ¼ β A ðβ A þ1Þ À ΓI β I þ1 ð Þ will always be positive, since Γ I 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.
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 (K P /K PI ) is small (≈0.08). 11 The Monod growth rate with inhibition is given by Eq. (1). Developing the equation further, it is found that When (K P /K PI ) ≪ 1, then Eq. (24) becomes 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, X f , and a constant average substrate effective diffusivity, D fP . In the present continuum heterogeneous biofilm model, the dimensionless effective diffusivity of phenol and biofilm density change with biofilm position (x Ã ) as: Therefore, the average biofilm density, X f , and the average substrate effective diffusivity, D fP , 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, C Ps , 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/m 3 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, C Ps < 50 g/m 3 , was due to the fact that the rate tends towards a global first order reaction with respect to phenol concentration. The inhibition effect ([K PI /(K PI + C Ps )] ≈ 1) was very small. However, when the phenol concentration was high (C Ps 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/m 3 , inhibition played an important role in the overall rate of phenol consumption. 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/m 3 ), 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 C Ps .
Finally, it is interesting to note that the effectiveness factor takes on values greater than one, for the three cases where C Ps > 200 g/ m 3 . To assure that η will take on values greater than one (case (c)) requires:   Another interesting observation regarding the behavior of η at different C Ps (Fig. 3b) is the following: for each biofilm thickness, the point where η = 1 corresponds to a C Ps 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 10 4 , 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 semianalytical 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/m 3 . 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 -C Pro /C ProCrit ), as suggested by Levenspiel,24 where C Pro is the product concentration and C ProCrit 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. 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, r ob , and b effectiveness factor, η, as a function of substrate biofilm-surface concentration, C PS , for different biofilm thicknesses, L f , of 75, 150 and 300 μm