Modeling Heavy Metal Sorption and Interaction in a Multispecies Biofilm

Abstract

A mathematical model able to simulate the physical, chemical and biological interactions prevailing in multispecies biofilms in the presence of a toxic heavy metal is presented. The free boundary value problem related to biofilm growth and evolution is governed by a nonlinear ordinary differential equation. The problem requires the integration of a system of nonlinear hyperbolic partial differential equations describing the biofilm components evolution, and a systems of semilinear parabolic partial differential equations accounting for substrates diffusion and reaction within the biofilm. In addition, a semilinear parabolic partial differential equation is introduced to describe heavy metal diffusion and sorption. The biosoption process modeling is completed by the definition and integration of other two systems of nonlinear hyperbolic partial differential equations describing the free and occupied binding sites evolution, respectively. Numerical simulations of the heterotrophic-autotrophic interaction occurring in biofilm reactors devoted to wastewater treatment are presented. The high biosorption ability of bacteria living in a mature biofilm is highlighted, as well as the toxicity effect of heavy metals on autotrophic bacteria, whose growth directly affects the nitrification performance of bioreactors.

1. Introduction

Most of the living microbial communities organize themselves in complex structures where the interaction between different species leads to advantageous environmental conditions for their growth [1,2]. These structures, known as multispecies biofilms, include different components, such as living cells, inert materials and extracellular polymeric substances (EPS) [3,4,5,6]. Their structural organization confers to these biological systems enhanced mechanical characteristics and adaptive features to many environmental conditions [7,8,9]. For instance, the protective self-secreted EPS matrix can strongly affect the dynamics of substances within the biofilm and it can also serve as a source of nutrients for bacteria [10,11].

These aspects are highly relevant in many applications as biofilms result in being more resistant than individual planktonic cells to toxic substances such as heavy metals, antibiotics, chlorine and detergents, due to the presence of natural diffusion barriers [12]. In recent years, biofilms have been widely used as biosorption technologies for metal immobilization and sequestration [13,14]. Biosorption is a combination of complex phenomena leading to the entrapment of a substance onto the surface of a living/dead organism or EPS. The mechanisms involved (complexation, precipitation, ion exchange, adsorption) are strongly affected by several biotic and abiotic parameters, such as pH, temperature, binding site density and affinity, which in turn influence the biosorption efficiency. Significant applications of biofilm technology to biosorption have been presented in the field of groundwater purification and mining industry wastewater treatment.

The sorption properties of various components constituting a biofilm (i.e., microorganisms, EPS and inert materials) depend on the different affinity of each specific component to heavy metals. It is known, for instance, that the cell membrane of many microorganisms allows for heavy metals accumulation due to the presence of surface functional groups [15]. These act as binding agents removing heavy metals during biofilm growth. On the other hand, heavy metals can be highly toxic compounds for a wide range of bacteria, i.e., autotrophic microorganisms, as they can act as inhibiting agent when significant metal concentrations are reached in bioreactors [16,17,18,19].

Many experimental studies demonstrated the possibility of using bacteria to govern heavy metal mobility in different aquatic ecosystems [20,21,22], but additional efforts are still required to completely understand the complex dynamics and interactions occurring between biofilms and heavy metals. In this context, mathematical modeling represents an appropriate tool to provide basic information on specific biosorption phenomena and stimulate further research on the multiplicity of mechanisms regulating biosorption process by biofilms [2]. For instance, multidimensional models can be implemented for specific applications when micro-scale outputs are required [7]. The spatial distribution of diffusing compounds and microbial species within the biofilm, and the physical structure of the biofilm at a micro-scale level can be investigated by using complex 2D and 3D mathematical models [23,24,25,26]. If a macro-scale output is required, as in the case of engineering biofilm reactors, 1D formulations have been recognized as efficient tools to analyze bioreactor performances in terms of biomass accumulation and degradation of substrates [27].

To this aim, a 1D mathematical model reproducing a biosorption phenomenon occurring in a typical biofilm reactor devoted to wastewater treatment has been proposed. The model is presented in its general form and then applied to a relevant case in wastewater treatment field. Specifically, the case study accounts for the coexistence of two different microbial species performing nitrification and organic carbon degradation. A continuum approach was used to describe biomass growth and decay within the biofilm [28]. The model accounts for the diffusion–reaction of substrates and the diffusion–biosorption of heavy metals within the biofilm [29,30]. More precisely, in this work, heterotrofic bacteria have been characterized by a high specific number of binding sites on their cell wall allowing heavy metal sequestration during biofilm evolution. On the other hand, the kinetics of autotrophic bacteria, which are usually more sensitive to toxic compounds then heterotrophic species, have been supposed to be negatively affected by the heavy metal concentration, which acts as an inhibiting agent and affects the efficiency of the nitrification process.

The main objective is to apply the knowledge of recently introduced mathematical approaches for biosorption in multispecies biofilms to highlight the effects of heavy metals in a traditional biofilm system for wastewater treatment. The work elucidates different ecological aspects of biofilms/heavy metals interaction, such as spatial distribution of biofilm components over time, substrate and heavy metals dynamics, and effects of heavy metals contamination. Numerical simulations remarked on the consistency of the model and showed the effect of toxic heavy metals on different microbial species coexisting in a multispecies biofilm.

2. Statement of the Problem

The effect of an inhibiting agent diffusing in a multispecies biofilm and the related biosorption interactions are discussed in the following sections. The specific case study concerns the competition for oxygen of heterotrophic and autotrophic microorganisms performing organic carbon degradation and nitrification, respectively. This is a typical situation occurring in the biological treatment units of municipal wastewater treatment plants.

According to [29], the biofilm dynamics were modeled as a free boundary problem essentially hyperbolic, where the free boundary is represented by the biofilm thickness. Its evolution is dictated by the growth of the microbial species constituting the biofilm $X_i(z,t)$ and the exchange fluxes between the biofilm and the bulk liquid. The biofilm growth is catalyzed by the availability of substrates $S_j(z,t)$, which diffuse from the bulk liquid into the biofilm where they are consumed by microbial species. The model considers biofilm as constituted of four different components $X_i,i=1,\dots ,4$ (green and gray in Figure 1), which can accumulate/growth and decrease/decay during time. The biofilm growth and development is governed by the availability of substrates $S_j,j=1,\dots ,3$ (blue in Figure 1) within the biofilm, which regulate the microbial metabolism and interactions. These components include heterotrophic bacteria $X_1={\rho }_1{f}_1$, autotrophic bacteria $X_2={\rho }_2{f}_2$, inert material $X_3={\rho }_3{f}_3$ and EPS $X_4={\rho }_4{f}_4$, with $f_i$ denoting the volume fraction of each biofilm component i and ${\rho }_i$ the corresponding density. Specifically, EPS production was taken into account according to the approach proposed by [31]. Three different substrates, ammonium $S_1$, organic carbon $S_2$, and oxygen $S_3$ were taken into account as they are involved in metabolic pathways. The active microbial biomasses $X_1$ and $X_2$ naturally decrease via respiration and decay processes, producing residual inert biomass $X_3$. Contextually, they produce EPS as a metabolic byproduct during their growth. The autotrophs $X_2$ are nitrifying bacteria that grow by consuming ammonium $S_1$ and oxygen $S_3$. On the other hand, the heterotrophic bacteria simultaneously uptake organic carbon $S_2$ and oxygen $S_3$ for their growth. The two species compete for space and oxygen in multispecies biofilms [28].

The inhibiting agent $\mu$ (orange in Figure 1) was assumed to interact with the biofilm in two different ways: it can adsorb on a specific biofilm component, e.g., the heterotrophic biomass $X_1$, and act as inhibiting agent for an active microbial biomass, e.g., the autotrophic bacteria $X_2$. Note that a single inhibiting agent $\mu$ was considered in this work, but, in more complex cases, the effect of different heavy metals $\mu$${}_k$, $k=1,\dots ,l$ can be taken into account by using a similar approach. The concentration of heavy metals in biofilm reactors negatively affects the kinetics of autotrophic bacteria, which are typically more sensitive to contamination than heterotrophic species. Consequently, a specific inhibition term was exclusively introduced in the autotrophic growth rate function. The sorption phenomenon was modeled by directly taking into account the dynamics of the binding sites of the biofilm matrix.

The biofilm growth is governed by the following equations:

$$\frac{\partial X_i}{\partial t}+\frac{\partial }{\partial z}\left(u{X}_i\right)={\rho }_i{r}_{M,i}^{}(z,t,\mu ,\mathbf{X},\mathbf{S}),i=1,\dots ,4,0\le z\le L\left(t\right),t>0,$$

$$X_i(z,0)=X_{i0}\left(z\right),i=1,\dots ,4,0\le z\le L_0,$$

(1)

$$\frac{\partial u}{\partial z}=\sum _{i=1}^4{r}_{M,i}(z,t,\mu ,\mathbf{X},\mathbf{S}),0<z\le L\left(t\right),t\ge 0,u(0,t)=0,t\ge 0,$$

(2)

$$\dot{L}\left(t\right)=u(L\left(t\right),t)+{\sigma }_a\left(t\right)-{\sigma }_d\left(L\left(t\right)\right),t>0,L\left(0\right)=L_0,$$

(3)

where $X_i={\rho }_i{f}_i(z,t)$ denotes the concentration of the four biofilm components considered, ${\rho }_i$ is the constant density, $u(z,t)$ is the velocity of microbial mass displacement with respect to the biofilm substratum, $r_{M,i}^{}(z,t,\mu ,\mathbf{X},\mathbf{S})$ is the biomass growth rate, $L\left(t\right)$ is the biofilm thickness, $\mathbf{X}=(X_1,X_2,X_3,X_4)$, and $\mathbf{S}=(S_1,S_2,S_3)$. Equation (1) is derived from local mass balance considerations and governs the growth of the microbial species constituting the biofilm. The biomass expansion is modelled as an advective flux and depends on the metabolic reactions carried out by the microbial species. The reaction terms $r_{M,i}$ account for the microbial growth and decay, and EPS production. Equation (2) governs the biomass growth velocity; it is obtained summing over i Equation (1) and considering the constrain ${\sum }_{i=1}^n{f}_i=1$. The biofilm thickness evolution is ruled by an ordinary differential equation (Equation (3)) that is derived from global mass balance considerations and depends on both the biomass growth velocity $u\left(L\right(t),t)$ and the detachment ${\sigma }_d\left(L\left(t\right)\right)$ and attachment ${\sigma }_a\left(t\right)$ fluxes [32,33]. The latter represent the exchange fluxes between the biofilm and the bulk liquid compartment.

The kinetic terms $r_{Mi}(z,t,\mu ,\mathbf{X},\mathbf{S})$ for the biofilm components $X_1$, $X_2$, $X_3$, and $X_4$ can be expressed as specified in the following lines. For the active biomass $X_1$ and $X_2$,

$$r_{M,1}=((1-k_1)K_{max,1}\frac{S_2}{K_{1,2}+S_2}\frac{S_3}{K_{1,3}+S_3}-b_{m,1}{F}_1\frac{S_3}{K_{1,3}+S_3}-(1-F_1)c_{m,1})X_1,$$

(4)

$$r_{M,2}=((1-k_2)K_{max,2}\frac{S_1}{K_{2,1}+S_1}\frac{S_3}{K_{2,3}+S_3}\frac{K_I}{K_I+\mu }-b_{m,2}{F}_2\frac{S_3}{K_{2,3}+S_3}-(1-F_2)c_{m,2})X_2,$$

(5)

and, for the inert component $X_3$,

$$r_{M,3}=(1-F_1)c_{m,1}{X}_1+(1-F_2)c_{m,2}{X}_2,$$

(6)

while, for the EPS component $X_4$,

$$r_{M,4}=k_1{K}_{max,1}\frac{S_2}{K_{1,2}+S_2}\frac{S_3}{K_{1,3}+S_3}{X}_1+k_2{K}_{max,2}\frac{S_1}{K_{2,1}+S_1}\frac{S_3}{K_{2,3}+S_3}\frac{K_I}{K_I+\mu }{X}_2,$$

(7)

where $K_{max,i}$ denotes the maximum growth rate for biomass i, $k_i^{}$ is the coefficient associated with EPS formation, $K_{i,j}$ represents the affinity constant of substrate j for biomass i, $b_{m,i}$ denotes the endogenous rate for biomass i, $c_{m,i}$ is the decay–inactivation rate for biomass i, $F_i^{}$ represents the biodegradable fraction of biomass i, $\mu$ is the concentration of the heavy metal, which is supposed toxic for autotrophic bacteria $X_2$, and $K_I$ is the inhibition constant. The kinetic growth rates for the inert material (Equation (6)) end EPS component (Equation (7)) are directly connected to the biological activities performed by the microbial species. These are modeled by Monod-like kinetics (Equations (4) and (5)) regulated by the availability of substrates within the biofilm.

The evolution of the free ${\vartheta }_i(z,t)$ and occupied $\overline{{\vartheta }_i}(z,t)$ binding sites is modeled by the equations

$$\frac{\partial {\vartheta }_i}{\partial t}+\frac{\partial }{\partial z}\left(u{\vartheta }_i\right)=r_{M,i}(z,t,\mu ,\mathbf{X},\mathbf{S})-r_{D,i}(z,t,\mu ,\mathbf{\vartheta }),i=1,\dots ,4,0\le z\le L\left(t\right),t>0,$$

$${\vartheta }_i(z,0)={\vartheta }_{i0}\left(z\right),i=1,\dots ,4,0\le z\le L_0,$$

(8)

$$\frac{\partial \overline{{\vartheta }_i}}{\partial t}+\frac{\partial }{\partial z}\left(u\overline{{\vartheta }_i}\right)=r_{D,i}(z,t,\mu ,\mathbf{\vartheta }),i=1,\dots ,4,,t>0,0\le z\le L\left(t\right),$$

$$\overline{{\vartheta }_i}(z,0)=\overline{{\vartheta }_{i0}}\left(z\right),i=1,\dots ,4,0\le z\le L_0,$$

(9)

where $r_{D,i}$ denotes the sorption rate, and ${\vartheta }_{i0}$ and $\overline{{\vartheta }_{i0}}$ are the initial distribution of the free and occupied binding sites, respectively. The free binding site fractions can increase (Equation (8)) due to the generation of new biomass, or decrease due to the biosorption. A parabolic partial differential equation (PDE) describes the evolution of the adsorbing compound $\mu$ within the biofilm

$$\frac{\partial \mu }{\partial t}-\frac{\partial }{\partial z}\left(D_{\mu }\frac{\partial \mu }{\partial z}\right)=-Y_{ADS}{N}_{\mu }{r}_D(z,t,\mu ,{\mathbf{\vartheta }}_i),0<z<L\left(t\right),t>0,$$

$$\mu (z,0)={\mu }_0\left(z\right),\frac{\partial \mu }{\partial z}(0,t)=0,\mu (L\left(t\right),t)={\mu }_L\left(t\right),0\le z\le L_0,t>0,$$

(10)

where $D_{\mu }$ is the diffusivity coefficient for the adsorbing compound, $N_{\mu }$ denotes the binding sites density and $Y_{ADS}$ is the yield of the adsorbing compound. The kinetic term $r_D$ describes a non-reversible heavy metal sorption mechanism. This is expressed by

$$r_D=K_{ads}\mu {\vartheta }_1,$$

(11)

where $K_{ads}$ denotes the adsorption constant. According to Equations (10) and (11), the dynamics of the adsorbing compound $\mu$ are regulated by the sorption rate $r_D$, which is multiplied by two parameters with physical meaning; $Y_{ADS}$ is the amount of adsorbing compound allocated in each binding site, and $N_{\mu }$ is the number of binding sites related to the specific biofilm component. These parameters describe the sequestration ability of a specific biofilm component.

The diffusion−reaction of each substrate was modeled by the equations

$$\frac{\partial S_j}{\partial t}-\frac{\partial }{\partial z}\left(D_{S,j}\frac{\partial S_j}{\partial z}\right)=r_{S,j}(z,t,\mu ,\mathbf{X},\mathbf{S}),j=1,\dots ,3,0<z<L\left(t\right),t>0,$$

$$S_j(z,0)=S_{j0}\left(z\right),\frac{\partial S_j}{\partial z}(0,t)=0,S_j(L\left(t\right),t)=S_{jL},j=1,\dots ,3,0\le z\le L_0,t>0,$$

(12)

where $D_{S,j}$ is the diffusivity coefficient, and $r_{S,j}(z,t,\mu ,\mathbf{X},\mathbf{S})$ is the conversion rate of substrate j. These terms are specifically expressed by

$$r_{{S,1}^{}}=-\frac{1}{Y_2}((1-k_2)K_{max,2}\frac{S_1}{K_{2,1}+S_1}\frac{S_3}{K_{2,3}+S_3}\frac{K_I}{K_I+\mu }{X}_2,$$

(13)

$$r_{{S,2}^{}}=-\frac{1}{Y_1}((1-k_1)K_{max,1}\frac{S_2}{K_{1,2}+S_2}\frac{S_3}{K_{1,3}+S_3}{X}_1,$$

(14)

$$r_{{S,3}^{}}=-(1-k_1)\frac{(1-Y_1)}{Y_1}((1-k_1)K_{max,1}\frac{S_2}{K_{1,2}+S_2}\frac{S_3}{K_{1,3}+S_3}{X}_1,$$

$$-(1-k_2)\frac{(1-Y_2)}{Y_2}((1-k_2)K_{max,2}\frac{S_1}{K_{2,1}+S_1}\frac{S_3}{K_{2,3}+S_3}\frac{K_I}{K_I+\mu }{X}_2,$$

$$-b_{m,1}{F}_1\frac{S_3}{K_{1,3}+S_3}{X}_1-b_{m,2}{F}_2\frac{S_3}{K_{2,3}+S_3}{X}_2,$$

(15)

where $Y_i$ denotes the yield for biomass i. A schematic representation of the model structure is shown in Figure 1.

3. Numerical Simulation

The presented mathematical model was applied to simulate the effect of exposition to a toxic heavy metal in a multispecies biofilm with an initial thickness of 300 $\mathsf{\mu }$m. The metal represents an adsorbing compound for one of the microbial species and acts as a toxic agent for the other. In particular, $\mu$ is supposed to be toxic for autotrophic bacteria but can be sorbed on the cellular membrane of heterotrophic bacteria. The values of the kinetic and stoichiometric parameters, and the mass transfer coefficients are reported in

Table 1. Kinetic parameters used for model simulations.

Parameter	Definition	Unit	Value
$K_{max,1}$	Maximum growth rate for $X_1$	d${}^{-1}$	4.8
$K_{max,2}$	Maximum growth rate for $X_2$	d${}^{-1}$	0.95
$k_1$	EPS formation by $X_1$	mg COD/mg COD	0.02
$k_2$	EPS formation by $X_2$	mg COD/mg COD	0.011
$K_{1,2}$	Organics half saturation constant for $X_1$	mg COD L${}^{-1}$	5
$K_{1,3}$	Oxygen half saturation constant for $X_1$	mg L${}^{-1}$	0.1
$K_{2,1}$	Ammonium half saturation constant for $X_2$	mg N L${}^{-1}$	1
$K_{2,3}$	Oxygen half saturation constant for $X_2$	mg L${}^{-1}$	0.1
$b_{m,1}$	Endogenous rate for $X_1$	d${}^{-1}$	0.025
$b_{m,2}$	Endogenous rate for $X_2$	d${}^{-1}$	0.0625
$F_1$	Biodegradable fraction of $X_1$	–	0.8
$F_2$	Biodegradable fraction of $X_2$	–	0.8
$c_{m,1}$	Decay-inactivation rate for $X_1$	d${}^{-1}$	0.05
$c_{m,2}$	Decay-inactivation rate for $X_2$	d${}^{-1}$	0.05
$Y_1$	Yield of $X_1$	g${}_{biomass}$/g${}_{substrate}$	0.4
$Y_2$	Yield of $X_2$	g${}_{biomass}$/g${}_{substrate}$	0.22
$Y_{ads}$	Yield of adsorbent	g${}_{metal}$/n${}_{sites}$	1
$N_{\mu }$	Binding sites density for $X_1$	n${}_{sites}$ L${}^{-1}$	1 and 50
$K_I$	Inhibition constant	mg L${}^{-1}$	${10}^{-5}$
$\rho$	Biofilm density	g m${}^{-3}$	2500
$\lambda$	Biomass shear constant	mm d${}^{-1}$	1250

. They were adopted according to [30]. The initial conditions and biological parameters adopted in the simulations are reported in

Table 2. Initial conditions for biofilm growth.

Parameter	Symbol	Unit	Value
COD concentration at $L=L\left(t\right)$	$S_{1L}$	mgL${}^{-1}$	20
Oxygen concentration at $L=L\left(t\right)$	$S_{3L}$	mgL${}^{-1}$	8
Ammonium concentration at $L=L\left(t\right)$	$S_{2L}$	mgL${}^{-1}$	2
Free metal concentration at $L=L\left(t\right)$	${\mu }_L$	mgL${}^{-1}$	2
Time Simulation	T	d	100
Initial Biofilm thickness	$L_0$	mm	0.3
Initial Volume Fraction of Autotrophs $\left(X_1\right)$	$f_{1,0}\left(z\right)$	–	0.399
Initial Volume Fraction of Heterotrophs $\left(X_2\right)$	$f_{2,0}\left(z\right)$	–	0.5
Initial Volume Fraction of Inert $\left(X_3\right)$	$f_{3,0}\left(z\right)$	–	0.001
Initial Volume Fraction of EPS $\left(X_4\right)$	$f_{4,0}\left(z\right)$	–	0.1

.

Numerical solutions to the free boundary problem stated in Section 2 were obtained by using the method of characteristics, e.g., [34,35,36,37]. Accuracy was checked by comparison to the geometric constraint ${\sum }_{i=1}^4{f}_i(z,t)=1$. Simulations were performed using an original software developed on the Matlab platform.

For all the dissolved species, i.e., substrates and adsorbing contaminant, Dirichlet conditions on the free boundary were assumed. In Equation (3) governing the free boundary evolution, ${\sigma }_d\left(L\left(t\right)\right)$ was assumed to be a known function of L and t

$${\sigma }_d\left(L\left(t\right)\right)=\lambda L_{}^2\left(t\right),$$

(16)

where $\lambda$ is the share constant whose value is reported in Table 1. No attachment phenomena were considered for all the numerical simulation, thus ${\sigma }_a\left(t\right)$ was fixed to zero. The initial biofilm composition is defined in Table 2. In particular, the biofilm is set to be initially constituted by the autotrophic (39.9%) and heterotrophic (50%) bacteria, EPS (10%) and inert (0.1%). The simulations reproduce a typical environmental condition occurring in the biological units of municipal wastewater treatment plants. The oxygen concentration in the bulk liquid was fixed to 8 mg/L, consistently with real scale continuous aerated systems. The concentrations of soluble organic carbon, i.e., chemical oxygen demand (COD), and ammonium, i.e., nitrogen content (N), in the bulk liquid were fixed to 20 mg/L and 2 mg/L, respectively.

The model outputs are reported in Figure 2 and Figure 3. Numerical simulations demonstrate model capability of predicting the spatial distribution of biofilm components, the occupied and free binding site fractions, the substrate trends, the free contaminants profiles over biofilm depth and the biofilm thickness. The simulations show the effect of the biosorption phenomenon on the biological evolution of the overall system, and how the different features of the heterotrophic biomass, such as the binding site density $N_{\mu }$, can substantially affect the final configuration of the biofilm and its properties.

Two different values of the binding site density $N_{\mu }$ were used for numerical simulations (Figure 2 and Figure 3). In the first simulation set, the value of site density $N_{\mu }$ was fixed to 1 (Figure 2), while, in the second set, $N_{\mu }$ was increased to 50 (Figure 3). In the first case ($N_{\mu }=1$), the low site density determines a low adsorption rate resulting in a high diffusion of the free contaminant, which shows a fully penetrated profile (Figure 2, C1). The presence of $\mu$ in the inner part of the biofilm inhibits the metabolic activities of autotrophic bacteria. Despite the presence of autotrophs into the biofilm (Figure 2, A1), the ammonium is not degraded; indeed, the concentration of ammonium in the system remains constant (Figure 2, B1). The fraction of autotrophic bacteria decreases with time due to the toxic effect of $\mu$ (Figure 2, A2, A3 and A4). After 100 days of simulation, the autotrophs completely disappear from the biofilm (Figure 2, A4).

In the second case ($N_{\mu }=50$), the higher site density determines a higher adsorption rate resulting in a lower diffusion of the heavy metal than in the first case (Figure 2, C1). It is interesting to notice that the concentration of the metal $\mu$ is essentially zero in the inner part of the biofilm, allowing the proliferation of autotrophic bacteria (Figure 3, A1). Due to the absence of $\mu$, the metabolic activity of autotrophic bacteria is not inhibited and the ammonium is degraded (Figure 3, B1). Notably, the existence of autotrophic bacteria in the inner part of the biofilm is due to the relevant adsorption phenomenon occurring in the external part of the biofilm. Indeed, heterotrophic bacteria act as a biological shield for autotrophic bacteria, which can live and proliferate in the biofilm structure performing their biological activity. The coexistence of the two species is preserved after 100 d of simulation as it is possible to notice from the final distribution of the microbial species within the biofilm (Figure 2, A4).

Additional simulations were run by varying the inhibition constant $K_I$ and the metal concentration within the bulk liquid ${\mu }_L$. The first set of simulations was performed to test the effect of an increasing resistance of the autrotrophic component $X_2$ to the toxic metal. The site density was set to $N_{\mu }=1$ to reproduce the case of a heterotrophic-autotrophic biofilm with a low sorption capability. The results were summarized in Figure 4.

When varying $K_I$ from ${10}^{-5}$ to ${10}^{-4}$ (Figure 4A), the autotrophic fraction rises from 0 to 15%. By further increasing the value of $K_I$ to ${10}^{-3}$ and ${10}^{-2}$, the autotrophic fraction reaches 18 and 17%, respectively. This slight difference is due to the biofilm thickness, which is smaller in the case of $K_I={10}^{-3}$, and affects the diffusion of substrates within the biofilm (Figure 4B). Ammonia shows a fully penetrated profile for all values of $K_I$, due to the low concentration of autotrophic bacteria (<20%) within the biofilm. Oxygen is characterized by a similar trend for the lowest values of $K_I$. When the autotrophic fraction increases as a result of a higher $K_I$ value, the oxygen concentration decreases all over the biofilm due to the additional consumption related to the autotrophic metabolism. No significant differences can be noted in the $\mu$ profile, which shows a fully penetrated profile for all $K_I$ values (Figure 4C). The simulation results highlight the key role played by the inhibition constant in the definition of the biofilm composition.

The second set of simulations was performed by varying the metal concentration in the bulk liquid to test the effectiveness of the biological shield provided by heterotrophic bacteria (Figure 5). The final simulation time and the site density were set at $T=100$ d and $N_{\mu }=50$, respectively. For ${\mu }_L$ = 4 ×${10}^{-2}$, the metal showed a fully penetrated profile (Figure 5C), which determines a strong inhibition of the autotrophic species. For all the other values of ${\mu }_L$, the metal concentration reaches zero within the biofilm.

Increasing ${\mu }_L$ from $4\times {10}^{-4}$ to $2\times {10}^{-2}$, the volume fraction of the autotrophic species slightly increases (Figure 5A) due to the small difference in biofilm thickness and substrates trends within the biofilm. When the autotrophic fraction is inhibited by the high metal concentration, ammonia remains constant within the biofilm and oxygen shows a fully penetrated profile. For all the simulations, the COD profile is invariant (Figure 5B). Except for ${\mu }_L=4\times {10}^{-2}$, the simulation results prove the effectiveness of the heterotrophic component in protecting the autotrophic fraction from metal exposure.

Further numerical simulations were carried out to test the influence of the initial distribution of biofilm components in both the experimental cases $N_{\mu }=1$ and $N_{\mu }=50$ (data not shown). After 100 days of simulation time, numerical results showed a negligible variation of the biofilm components distribution and a similar biological response to the heavy metal exposition.

4. Conclusions

In this work, a free boundary problem related to biofim growth and evolution during heavy metal exposition in wastewater treatment plants has been discussed. The model highlights the dynamic interactions occurring between different biofilm components when an inhibiting compound diffuses from the bulk liquid within the biofilm structure. The biosorption phenomenon has been considered by assigning a specific binding site density to the heterotrophic biomass, which is able to act as a biological shield for autotrophic bacteria. Numerical results showed the crucial role of heterotrophic bacteria on biosorption processes occurring in wastewater treatment plants. The combined effect of heavy metals inhibition and biosorption phenomena within the biofilm structure has been newly analyzed in this study. The general form and the structure of the mathematical model allow for its application to different biological cases of high engineering interest. Simulation results demonstrated that biofilm systems can be effectively used in the context of bioremediation. The development of 1D mathematical models able to predict biofilm evolution and features under different environmental condition is highly relevant for real scale applications, such as heavy metal recovery in biofilm reactors. Further experimental studies are still required to elucidate the different interactions occurring between heavy metals and specific biofilm components. For instance, the role of EPS on heavy metals biosorption can affect the biological response of a specific multispecies biofilm. This role could be further taken into account with the presented mathematical model by assigning a specific site density for the EPS component as a function of the biosorption affinity with the diffusing heavy metal.