Modeling of Biofilm Growth on Fine Spherical Particles with the Use of Cellular Automata: The Influence of Cell Death and Lysis on the Biofilm Structure

The paper concerns the modeling of heterogeneous biofilm growth on fine spherical particles of such biofilm forms as, e.g., fluidized-bed bioreactors. Three discrete mathematical models based on cellular automata theory were proposed. The double-substrate kinetics of biomass growth, biomass displacement, internal and external mass transfer resistances, death and lysis of microbiological cells and biofilm detachment were taken into account. It was shown that there are no significant qualitative and quantitative differences between biofilm growth on flat and spherical particles of different radii. Computer simulations were compared with experimental observations. Qualitative and quantitative agreement areachieved if both cell death and lysis aretaken into consideration and a proper algorithm of biomass displacement is used. The value of the bacteria lysis rate coefficient was estimated.


Introduction
The biological processes of wastewater treatment can be divided into suspended-, attachedand mixed-growth ones. In suspended-growth processes, microorganisms flow freely in the liquid phase. Microorganisms in attached-growth processes grow on solid substratum in theform of a biofilm, whereas, in mixed-growth processes, both forms are present in a system. Particle-based biofilm reactors, in which fine particles are used as a substratum for biofilm growth, have been widely employed in numerous processes, including the biodegradation of pollutants [1] and the production of antibiotics [2]. Reactors belonging to this group; for example, fluidized-bed bioreactors and air-lift suspension bioreactors can be operated at high volumetric flow rates, since the biomass is immobilized on a solid carrier of highly developed interphase surface, and retained inside the reactor [3].
Numerous studies have been published during last five decades concerning modeling and understanding the processes occurring in biofilms. The first mathematical models of biofilms assumed aflat geometry and a uniform distribution of biomass [4]. Those models were used for describing the utilization and diffusion of a single limiting substrate. Advances in experimental techniques, such as microelectrode technology and confocal scanning laser microscopy (CSLM), enabled researchers to determine the structure of biofilm more accurately and to verify the results of computer simulations [5]. Discoveries concerning heterogeneous biofilm growth induced an interest in the modeling of biofilm structure. The effect of studies in this area are two-and three-dimensional mathematical models based on a properly defined grid (cellular automata) [6][7][8] andon individuals (individual-based models) [9].
Several cellular automata-based models of the biofilm growth have been published so far. However, no model is able to predict the morphology and structure of biofilm growing on fine spherical particles.
One of the few studies in this area is the work by Graf von Schulenburg et al. [10], who proposed a discrete model of biofilm growth in porous media. However, the quoted authors accepted the single-substrate Monod kinetic model and neglected biofilm detachment.
In this study, three models based on cellular automata are proposed. The first one neglects death and the lysis of microbiological cells and a previously published algorithm of biomass displacement is used. The second one takes into account the death of microorganisms, while the third one allows for the death and lysis of microorganisms. Moreover, a new algorithm of biomass displacement is proposed. Dynamical simulations of biofilm growth were performed using the proposed models and compared. The relationships between average biofilm density and its thickness were compared with those obtained experimentally.

Mathematical Model of the Biofilm Dynamics
The mathematical model is based on the two-dimensional cellularautomata model of biofilm growth on flat substratum presented in a previous study [11,12]. It was modified in order to represent growth on spherical particles. The model takes into account the diffusion and utilization of two growth-limiting substrates, the spreading of biomass and biofilm detachment. The nextmodification was introduced in order to take into account death and the lysis of bacterial cells.
The cellular automata-based models differ in their approach tosimulating biomass spreading and cell death and lysis. Appendix A presents rules which were used in all models, namely"diffusion of the substrates", "utilization of the substrates in the biofilm", "growth of microorganisms and biofilm detachment".
The formulated cellular automata differ in terms of the algorithm of biomass spreading used in the following ways: 1.
The algorithm of biomass spreading proposed by Picioreanu et al. [7]. Cell death and lysis are neglected. It will be further referred to as CA-1.

2.
The algorithm of biomass spreading proposed in this study. Cell death is taken into account. It will be further referred to as CA-2.

3.
An algorithm of biomass spreading that is the same as in CA-2. However, both death and lysis of microbiological cells are taken into account. It will be further referred to as CA-3.
The CA-1 model consists of three grids describing the concentration of biomass, carbonaceous substrate and oxygen. Due to the negligence of microbiological cell death, all biomass is active. Cellular automata CA-2 and CA-3 consist of four grids describing, respectively, concentrations of active biomass, dead (inert) biomass, carbonaceous substrate and oxygen. The cell states of each grid of the cellular automata correspond to ρ ba , ρ bo , c b A and c b T , respectively. In order to explain the algorithm for the death of microorganisms, let us consider the following kinetic expression of the biomass death rate: At every time step ∆t g , the biomass density will decrease due to the death of bacteria: where i and j are indices of the grid.
The death of active bacteria is connected with the generation of dead (inactive) microbiological cells. Therefore, the cell state of the grid representing the concentration of dead bacteria will change according to: Experimental studies indicate that inert biomass loses its mechanical properties [13]. The process of bacterial lysis is modeled in this study as a first-order process. Only dead cells are subject to this phenomenon: A similar approach was used in the cellular automata-based model proposed by Pizarro et al. [14]; however, this model used single-substrate kinetics and flat biofilm.
Algorithm for Biomass Spreading Used in CA-2 and CA-3 According to Tang and Valochi [15] the application of the algorithm proposed by Picioreanu et al. [7] results in significant discontinuities in the biomass distribution. Instead of a probabilistic algorithm for the biomass spreading, Tang and Valochi proposed a deterministic one, based on the determination of theshortest path from the source grid cell to the destination cell. Based on the observations of the quoted authors, the following algorithm was formulated.

1.
Checkif the following condition is fulfilled: If yes, go to step 2. Otherwise, finish the algorithm.

2.
Calculate the amount of biomass thatexceeds the maximum density: 3.
Calculate the amount of active and dead bacteria in the value calculated in point 1, as follows:

4.
Determine the concentrations of active and dead bacteria in time t + ∆t:

Dynamics of the Biofilm Growth
The computer program implementing the CA-1 model generates results in the form of three two-dimensional arrays at each time moment. Each array represents concentrations of biomass and two liming substrates: carbonaceous substrate A and oxygen T. For mathematical models thattake into account dead cells, i.e., CA-2 and CA-3, an additional array is created to represent dead bacteria. The value of each cell in a grid, i.e., the cell state, indicates the concentration of a given component.
The aerobic degradation of phenol was chosen as a process example due to the large significance of this process in wastewater treatment. Moreover, several experimental studies can be found in the literature, which can be used for the validation of model predictions [16][17][18]. The kinetic parameters of this process are presented in Table 1. Table 1. Kinetic parameters of aerobic degradation of phenol by Pseudomonas putida [19].  Figure 1 were obtained with the use of CA-1, i.e., with the negligence of cell death and lysis. It was accepted that, at an initial time, the particle is uniformly covered with biomass of density equal to 0.1 · ρ max b . The concentrations of phenol and oxygen in the bulk liquid were equal to 0.02 kg/m 3 and 0.008 kg/m 3 , respectively. The chosen dissolved oxygen concentration is close to the saturation state.

Dynamics of the Biofilm Growth
The computer program implementing the CA-1 model generates results in the form of three twodimensional arrays at each time moment. Each array represents concentrations of biomass and two liming substrates: carbonaceous substrate A and oxygen T. For mathematical models thattake into account dead cells, i.e., CA-2 and CA-3, an additional array is created to represent dead bacteria. The value of each cell in a grid, i.e., the cell state, indicates the concentration of a given component.
The aerobic degradation of phenol was chosen as a process example due to the large significance of this process in wastewater treatment. Moreover, several experimental studies can be found in the literature, which can be used for the validation of model predictions [16][17][18]. The kinetic parameters of this process are presented in Table 1. Table 1. Kinetic parameters of aerobic degradation of phenol by Pseudomonas putida [19].  Figure 1 were obtained with the use of CA-1, i.e., with the negligence of cell death and lysis. It was accepted that, at an initial time, the particle is uniformly covered with biomass of density equal to  It can be seen in Figure 1 that oxygen was depleted in deeper parts of the biofilm (Figure 1i-l), while phenol penetrated the whole biofilm (Figure 1e-h). Such behavior is related to the large requirement foroxygen in this process, which is reflected by yield coefficients wBA and wBT. According to these values, the utilization of 1 kg of phenol is related to the utilization of 1.54 kg of oxygen. Even if the biofilm thickness is relatively low, approximately 100 μm (Figure 1c), the biofilm is not fully penetrated by oxygen.
The biomass distributions presented in Figure 1 show that biofilm thickness may decrease during the process. Thisis due to sloughing, i.e., a phenomenon based on the sudden detachment of large portions of the biofilm. It occurs when external parts of the biofilm lose their connection with deeper layers attached to the substratum. The sloughing event is illustrated more clearly in Figure 2, where the relationship between maximum biofilm thickness and time is shown. This mechanism occurs twice during 200 h of the process. At a growth time equal to 86 h, the biofilm thickness decreased from approximately 121 μm to 87 μm, while, at 145 h, sloughing occurred for the second time.  It can be seen in Figure 1 that oxygen was depleted in deeper parts of the biofilm (Figure 1i-l), while phenol penetrated the whole biofilm (Figure 1e-h). Such behavior is related to the large requirement foroxygen in this process, which is reflected by yield coefficients w BA and w BT . According to these values, the utilization of 1 kg of phenol is related to the utilization of 1.54 kg of oxygen. Even if the biofilm thickness is relatively low, approximately 100 µm (Figure 1c), the biofilm is not fully penetrated by oxygen.
The biomass distributions presented in Figure 1 show that biofilm thickness may decrease during the process. Thisis due to sloughing, i.e., a phenomenon based on the sudden detachment of large portions of the biofilm. It occurs when external parts of the biofilm lose their connection with deeper layers attached to the substratum. The sloughing event is illustrated more clearly in Figure 2, where the relationship between maximum biofilm thickness and time is shown. This mechanism occurs twice during 200 h of the process. At a growth time equal to 86 h, the biofilm thickness decreased from approximately 121 µm to 87 µm, while, at 145 h, sloughing occurred for the second time. It can be seen in Figure 1 that oxygen was depleted in deeper parts of the biofilm (Figure 1i-l), while phenol penetrated the whole biofilm (Figure 1e-h). Such behavior is related to the large requirement foroxygen in this process, which is reflected by yield coefficients wBA and wBT. According to these values, the utilization of 1 kg of phenol is related to the utilization of 1.54 kg of oxygen. Even if the biofilm thickness is relatively low, approximately 100 μm (Figure 1c), the biofilm is not fully penetrated by oxygen.
The biomass distributions presented in Figure 1 show that biofilm thickness may decrease during the process. Thisis due to sloughing, i.e., a phenomenon based on the sudden detachment of large portions of the biofilm. It occurs when external parts of the biofilm lose their connection with deeper layers attached to the substratum. The sloughing event is illustrated more clearly in Figure 2, where the relationship between maximum biofilm thickness and time is shown. This mechanism occurs twice during 200 h of the process. At a growth time equal to 86 h, the biofilm thickness decreased from approximately 121 μm to 87 μm, while, at 145 h, sloughing occurred for the second time.  The occurrence of sloughing is related to the existence of voids in the middle of the biofilm (Figure 1b-d). This property is caused by the mutual influence of the biofilm detachment (shear) and oxygen depletion. According to the accepted detachment model [20], the probability of this phenomenon is proportional to the square of the distance from the carrier surface. The biomass, which is placed near the biofilm surface, is more susceptible to the detachment;however, the substrates are available in larger amounts, which facilitates biomass growth. In turn, the growth of biomass near the carrier substratum is limited due to the depletion of oxygen, but the detachment probability is low due to the small distance from the carrier surface.
The radii of sand particles used in fluidized-bed bioreactors vary in the range of about 100-400 µm [3]. Other materials can also be used as a substratum for biofilm growth. For example, Seker et al. [21] used carbon particles with a radius equal to 250 µm. It is justified to perform an analysis of the influence of the size of an inert particle on the biofilm structure. Computer simulations were performed for spherical particles with radiibetween 150 µm and 350 µm and for a flat substrata, for comparison. The relationships between the biofilm thickness and time for all simulated biofilms arepresented in Figure 2. Close qualitative and quantitative agreement can be seen. In every case, the sloughing occurs after about 80 h at abiofilm thickness of approximately 125 µm. The phenomenon occurs for the second time when the biofilm thickness reaches approximately 125 µm once again. Such fluctuations inthe biofilm thickness have been observed in laboratory experiments, for example, in the study by Horn et al. [22]. The phenomenon was called a "quasi-steady state" by the quoted authors. Figure 3 presents two-dimensional distributions of active and inert biomass obtained with the use of CA-2 for time moments equal to t = 50 h (Figure 3a,e, i.e., first row), t = 100 h (Figure 3b,f, i.e., second row), t = 150 h (Figure 3c,g, i.e., third row) and t = 200 h (Figure 3d,h, i.e., fourth row). The value of thecell death rate coefficient was accepted according to Kumar et al. [23], i.e., k o = 5.6 × 10 −3 1/h. It can be seen in Figure 3a,e that at t = 50 h, despite relatively large biofilm thickness, almost all bacteria are active, i.e., their concentrations are close to ρ max b . The dead bacteria concentrations are significantly lower. Thisis caused by the large difference between the biomass growth rate and death rate, i.e.,the value of k is approximately 100 times larger than k o . The biofilm growth rate is relatively large; therefore, 50 h of growth is sufficient for obtaining athickness over 80 µm. Due to the significantly lower death rate coefficient, inactive bacteria did not accumulate due to the short growth time. Inactive bacteria accumulates near the substratum, i.e., at the biofilm bottom. This is related to the utilization of substrates during transfer into the biofilm. Bacteria located at the carrier surface are subject to the lowest substrate concentrations; therefore, the growth is limited in this region. The distributions of inactive bacteria after larger values of time (Figure 3f-h) show that concentrations of inactive bacteria continuously increase at the biofilm bottom. Aggregation in this region is also related to a low detachment probability near the carrier surface. In turn, dead bacteria near the biofilm external surface are continuously removed due to the detachment. The depletion of substrates in deeper biofilm partsalso results in lower active biomass concentrations. The occurrence of sloughing is related to the existence of voids in the middle of the biofilm (Figure 1b-d). This property is caused by the mutual influence of the biofilm detachment (shear) and oxygen depletion. According to the accepted detachment model [20], the probability of this phenomenon is proportional to the square of the distance from the carrier surface. The biomass, which is placed near the biofilm surface, is more susceptible to the detachment;however, the substrates are available in larger amounts, which facilitates biomass growth. In turn, the growth of biomass near the carrier substratum is limited due to the depletion of oxygen, but the detachment probability is low due to the small distance from the carrier surface.
The radii of sand particles used in fluidized-bed bioreactors vary in the range of about 100-400 μm [3]. Other materials can also be used as a substratum for biofilm growth. For example, Seker et al. [21] used carbon particles with a radius equal to 250 μm. It is justified to perform an analysis of the influence of the size of an inert particle on the biofilm structure. Computer simulations were performed for spherical particles with radiibetween 150 μm and 350 μm and for a flat substrata, for comparison. The relationships between the biofilm thickness and time for all simulated biofilms arepresented in Figure 2. Close qualitative and quantitative agreement can be seen. In every case, the sloughing occurs after about 80 h at abiofilm thickness of approximately 125 μm. The phenomenon occurs for the second time when the biofilm thickness reaches approximately 125 μm once again. Such fluctuations inthe biofilm thickness have been observed in laboratory experiments, for example, in the study by Horn et al. [22]. The phenomenon was called a "quasi-steady state"by the quoted authors. Figure 3 presents two-dimensional distributions of active and inert biomass obtained with the use of CA-2 for time moments equal to t = 50 h (Figure 3a,e, i.e., first row), t = 100 h (Figure 3b,f, i.e., second row), t = 150 h (Figure 3c,g, i.e., third row) and t = 200 h (Figure 3d,h, i.e., fourth row). The value of thecell death rate coefficient was accepted according to Kumar et al. [23], i.e.,ko = 5.6•× 10 −3 1/h. It can be seen in Figure 3a,e that at t = 50h, despite relatively large biofilm thickness, almost all bacteria are active, i.e., their concentrations are close to max b ρ . The dead bacteria concentrations are significantly lower. Thisis caused by the large difference between the biomass growth rate and death rate, i.e.,the value of k is approximately 100 times larger than ko. The biofilm growth rate is relatively large; therefore, 50 h of growth is sufficient for obtaining athickness over 80 μm. Due to the significantly lower death rate coefficient, inactive bacteria did not accumulate due to the short growth time. Inactive bacteria accumulates near the substratum, i.e., at the biofilm bottom. This is related to the utilization of substrates during transfer into the biofilm. Bacteria located at the carrier surface are subject to the lowest substrate concentrations;therefore, the growth is limited in this region. The distributions of inactive bacteria after larger values of time (Figure 3f-h) show that concentrations of inactive bacteria continuously increase at the biofilm bottom. Aggregation in this region is also related to a low detachment probability near the carrier surface. In turn, dead bacteria near the biofilm external surface are continuously removed due to the detachment. The depletion of substrates in deeper biofilm partsalso results in lower active biomass concentrations.  Similar simulations were performed using the CA-3 model, i.e.,bytaking into account the death and lysis of microbiological cells. Due to a lack of experimental data concerning this phenomenon, the value of klys = 1.12 h −1 was used, i.e., twice larger than the death rate coefficient. Simulations later in the study were performed for different values of this coefficient. Two-dimensional distributions of active and inert biomass are presented in Figure 4 in a similar manner as in Figure 3. At the initial stage of biofilm formation, bothdistributions of active biomass (Figure 4a and Figure 4e) are very similar to thosepresented in Figure 3a,e, respectively. However, for longer growth times, significantly lower inert biomass concentrations are obtained when lysis is taken into consideration. Thiscan be seen in Figure 3fh. Surprisingly, thisdoes not affect low active biomass concentrations at the biofilm bottom. Thiscan be explained by the depletion of oxygen, and therefore the lack of biomass growth in this biofilm region. Similar simulations were performed using the CA-3 model, i.e., bytaking into account the death and lysis of microbiological cells. Due to a lack of experimental data concerning this phenomenon, the value of k lys = 1.12 h −1 was used, i.e., twice larger than the death rate coefficient. Simulations later in the study were performed for different values of this coefficient. Two-dimensional distributions of active and inert biomass are presented in Figure 4 in a similar manner as in Figure 3. At the initial stage of biofilm formation, bothdistributions of active biomass (Figure 4a,e) are very similar to thosepresented in Figure 3a,e, respectively. However, for longer growth times, significantly lower inert biomass concentrations are obtained when lysis is taken into consideration. Thiscan be seen in Figure 3f-h. Surprisingly, thisdoes not affect low active biomass concentrations at the biofilm bottom. Thiscan be explained by the depletion of oxygen, and therefore the lack of biomass growth in this biofilm region.

Comparison of Model Predictions with Experimental Observations
The results of several experiments of phenol degradation by Pseudomonas putida conducted in fluidized-bed bioreactors can be found in the literature. According to the results of these studies, the relationship 1/ρ b = f (L b ) is a straight line, andis presented in Figure 5. Thusfar, no model has been published which predicts such biomass distribution. It should be noted that the results of the quoted experiments were conducted in different process conditions, whichexplains the different slopes of the obtained lines. The density distribution in the biofilm depends on, among other things, the type of fluidized bed (i.e., two-phase or three-phase), the number of particles in the bed and the species of bacteria [24]. The presence of bubbles in the bed results in higher detachment rates. This phenomenon is also intensified by larger particle numbers due to an increased frequency of collisions.

Comparison of Model Predictions with Experimental Observations
The results of several experiments of phenol degradation by Pseudomonas putida conducted in fluidized-bed bioreactors can be found in the literature. According to the results of these studies, the is a straight line, andis presented in Figure 5. Thusfar, no model has been published which predicts such biomass distribution. It should be noted that the results of the quoted experiments were conducted in different process conditions, whichexplains the different slopes of the obtained lines. The density distribution in the biofilm depends on, among other things, the type of fluidized bed (i.e., two-phase or three-phase), the number of particles in the bed and the species of bacteria [24]. The presence of bubbles in the bed results in higher detachment rates. This phenomenon is also intensified by larger particle numbers due to an increased frequency of collisions.  Table 2. It can be seen that both the cell death rate coefficient and cell lysis rate coefficient influence the linearity of the The highest value of R 2 was obtained for klys = 1.12•× 10 −2 1/h. The obtained function is presented in Figure 6. It canalso be seen in Table 2 that the biomass spreading algorithm proposed by Picioreanu et al. [7] is characterized by the lowest R 2 value. According to the authors of this study, the mentioned biomass spreading algorithm should not be used for the modeling of biofilm growth. The results obtained confirm the conclusions drawn by Tang and Valocchi [15] that the probabilistic biomass spreading algorithm results in unrealistic biomass distributions.  Table 2. It can be seen that both the cell death rate coefficient and cell lysis rate coefficient influence the linearity of the relationship 1/ρ b = f (L b ). The highest value of R 2 was obtained for k lys = 1.12 × 10 −2 1/h. The obtained function is presented in Figure 6. It canalso be seen in Table 2 that the biomass spreading algorithm proposed by Picioreanu et al. [7] is characterized by the lowest R 2 value. According to the authors of this study, the mentioned biomass spreading algorithm should not be used for the modeling of biofilm growth. The results obtained confirm the conclusions drawn by Tang and Valocchi [15] that the probabilistic biomass spreading algorithm results in unrealistic biomass distributions.

Conclusions
Three cellular automata-based models of biofilm growth are proposed in thisstudy. The first one neglects the death and lysis of microbiological cells, and uses a previously published algorithm of biomass displacement. The second one takes into account the death of microorganisms, while the third one allows for the death and lysis of microorganisms. A new algorithm of biomass displacement is also proposed.
The obtained results prove that both the death and lysis of microorganisms significantly influence the biofilm structure. The models were validated by a comparison of the relationships between average biomass density and biofilm thickness, obtained numerically, with those obtained experimentally by Seker et al. [19]. The results of phenol degradation by Pseudomonas putida in fluidizedbed bioreactors were used for this purpose. It was shown that, in order to obtain relationships corresponding to the experimental ones, it is necessary to use a proper approach to biomass displacement and to take into account the death and lysis of microorganisms. Moreover, the value of the lysis rate coefficient was estimated.

Conclusions
Three cellular automata-based models of biofilm growth are proposed in thisstudy. The first one neglects the death and lysis of microbiological cells, and uses a previously published algorithm of biomass displacement. The second one takes into account the death of microorganisms, while the third one allows for the death and lysis of microorganisms. A new algorithm of biomass displacement is also proposed.
The obtained results prove that both the death and lysis of microorganisms significantly influence the biofilm structure. The models were validated by a comparison of the relationships between average biomass density and biofilm thickness, obtained numerically, with those obtained experimentally by Seker et al. [19]. The results of phenol degradation by Pseudomonas putida in fluidized-bed bioreactors were used for this purpose. It was shown that, in order to obtain relationships corresponding to the experimental ones, it is necessary to use a proper approach to biomass displacement and to take into account the death and lysis of microorganisms. Moreover, the value of the lysis rate coefficient was estimated.

Rule 1-Diffusion of the Substrates
The probabilities of the movement of the mass quantum of oxygen in the liquid phase in all four directions on the grid are the same and equal to P dTw , and the sum of these probabilities is equal to one. The probability of the movement of the mass quantum of oxygen is dependent on the presence of biomass in the starting and the target cell. It was calculated as follows [12]: The above equation is useful for both the water and the biofilm. In Equation (A1), n b is the number characterizing the system of two cells: the starting and the target one. Thisspecifies the sum of the cells in this doublet, which represent the biomass. In turn, n w is the number related to the presence of water in the same doublet of cells. Both n b and n w can take values from the set {0,1,2}.
Similarly, as was done above for the oxygen transfer process, the relationship defining the probability P dA of the displacement of the mass quantum of the carbonaceous substrate is as follows: It was assumed that the ratio of the probability of the mass quantum displacement of the carbonaceous substrate in water to the probability of the mass quantum displacement of oxygen in the water is equal to the ratio of diffusion coefficients in water, P dAw /P dTw = D Aw /D Tw .

Rule 2-Utilization of the Substrates in the Biofilm
The rates of utilization of the carbonaceous substrate r A and oxygen r T in the biofilm are described by Equations (A3) and (A4) The forms of functions f 1 and f 2 depend on a given microbiological process. For the aerobic degradation of phenol, they take the form of Equations (A5) and (A6), respectively: Using the equation defining the reaction rate related to substrate A, we have Based on the above relationship, the expressions determining the increases in the concentrations of substrates during the time step ∆t can be obtained. These expressions are shown below: A (k, l, t) = −∆t · r b A c b A (k, l, t), c b T (k, l, t), ρ ba (k, l, t) (A8) ∆c b T (k, l, t) = −∆t · r b T c b A (k, l, t), c b T (k, l, t), ρ ba (k, l, t) (A9) The new values of the cells' states in the grids for the substrates are c b,j A (k, l, t + ∆t) = c b,j A (k, l, t) + ∆c b A (k, l, t) (A10) c b,j T (k, l, t + ∆t) = c b,j T (k, l, t) + ∆c b T (k, l, t) (A11)

Rule 3-Growth of Microorganisms and Biofilm Detachment
Due to the significantly different values of the characteristic times of diffusion and utilization and those of the growth of biomass and biofilm detachment [9], this rule is performed for larger time step valuesin theorder of hours [11].
The increase in the mass of active microorganisms in one time step ∆t g is determined using Equation (A12) ∆ρ ba (k, l, t) = ∆t g · r b B c b A (k, l, t), c b T (k, l, t), ρ ba (k, l, t) (A12) The new value of the cell's state is then ρ ba k, l, t + ∆t g = ρ ba (k, l, t) + ∆ρ ba (k, l, t) To calculate the probability of the biofilm detachment, Equation (A14) was used, as proposed by Chambless et al. [20].