Interspecific competition for substrate and space gives rise to considerable variation in biomass distribution within the microbial community. To study microbial community in depth, we used several research methods as sampling and analytical measurements, and developed a cellular automata (CA) model that would facilitate a description of the microbial growth process based on Anaerobic Digestion Model No. 1 (ADM1) of the International Water Association (IWA). Using the CA model, we aimed to determine whether interspecific competition occurs among acidogens, acetogens and methanogens, and to examine the influence of interspecific competition on the spatial structure of microbial communities. We found that acetogens and methanogens competed for core space, resulting in a multi-layer structure. Butyrate-degrading acetogens increased in number, resulting in inhibition of propionate-degrading acetogens. Hydrogenotrophic methanogens showed stronger competitive advantage than acetotrophic methanogens. The simulation showed that the multi-layer structure of the microbial community was formed by interspecific competition.

  • A two-dimensional model is developed based on cellular automata and ADM1.

  • New biomass-spreading rules are introduced to cellular automata.

  • The integrated model is evaluated by simulation of ABR start-up.

  • The new model can be used to simulate interspecific competition.

Graphical Abstract

Graphical Abstract
Graphical Abstract

Anaerobic digestion has become a favored method of wastewater treatment because of its low production of biomass, high cost–effectiveness and high energy production (Ahmed & Rodríguez 2017; Cook et al. 2017). The formation of anaerobic granular sludge in anaerobic digesters allows considerably higher loading rates than those typically applied in conventional activated sludge processes (Gagliano et al. 2017; Baeten et al. 2019). The resulting reduction in the reactor size and required area for treatment leads to a lower investment cost (Hulshoff Pol et al. 2004). Sludge granules were first observed with the anaerobic filter by Young & McCarty (1969). Now, over 50 years later, numerous researchers from all over the world have studied the granulation process (Batstone et al. 2004; Picioreanu et al. 2005; Odriozola et al. 2016; Feldman et al. 2017). However, anaerobic granulation technology is still not used at its full potential because of the high complexity of the granulation process and lack of knowledge about the underlying mechanism (Grotenhuis et al. 1991; Fang et al. 1995). Over the past two decades tremendous research progress has been made in understanding the microbiological characteristics of granules and the interaction among different microbial species within the granules (Fang 2000). Most theories have indicated that the methanogens Methanosaeta and Methanobacterium strain AZ play a key role in granular growth (Sanchez et al. 1994; Wu et al. 1996). Microcolonies consisting of syntrophic propionate, or butyrate degraders and hydrogen-utilizing methanogens have been observed in granular sludge (Schmidt & Ahring 1996). In addition to the observed microcolonies in the granules (where the degradation of propionate or butyrate may occur), thermodynamic and flux considerations have shown that the most effective degradation of propionate and butyrate occurs in microcolonies in which the distance between the syntrophic bacteria is small (Pauss et al. 1990). Based on microscope observations, a multi-layer model of microbial community was initially proposed by MacLeod et al. (1990) and Guiot et al. (1992). Microbial distribution in granules strongly depends on the degradation thermodynamics and kinetics of each individual substrate. It appears that interspecific competition for substrate and space may give rise to different structural granules.

Cellular automata (CA) modeling can be applied to describe the formation of microcolonies (Laspidou & Rittmann 2004). The CA model is defined as a spatially and temporally discrete system in which the state of an automaton is determined by a set of rules that act locally but apply globally (Picioreanu et al. 1998). This approach avoids the mathematical and computational burden of continuous models. (Tang & Valocchi 2013). The CA model can produce a large variety of distinct morphologies in response to changes in growth conditions (Manukyan et al. 2017). Recently, based on the CA theory, a series of multi-dimensional models with heterogeneous biomass and substrate distribution in two or three dimensions has been developed (Laspidou et al. 2010). The Anaerobic Digestion Model No. 1 (ADM1), developed by the IWA Task Group for mathematical modeling of anaerobic digestion processes, is one of the most sophisticated anaerobic digestion models, involving 19 biochemical processes and two physiochemical processes (Batstone et al. 2002; Barrera et al. 2015). The hydrolysis, acidogenesis and methanogenesis are considered in ADM1. The ADM1 model is included in the standard model libraries in most software packages developed for simulation of wastewater treatment plants (WWTPs) (GPS-X, Mike-WEST and Simba) (Jeppsson et al. 2016). The ADM1 was also designed to be readily extendible from one dimension to multiple dimensions for modeling the microbial community (Batstone et al. 2006; 2015). The aims of this study were to: (i) describe the microbial growth process during anaerobic digestion based on an integrated model combining CA and ADM1; and (ii) evaluate the influence of interspecific competition on the microbial community structure.

Reactor and its operation

The four-compartment (sequentially named C1, C2, C3 and C4) anaerobic baffled reactor (ABR) was constructed by 10 mm thick transparent Perspex with size of 74.0, 56.0, and 11.0 cm in length width and height, respectively and an effective volume of 28.4 L (Figure 1). Each of the four identical compartments included an up-flow zone and a down-flow zone with a volume ratio of 4.6:1. A 1.5 L water tank was attached to the last compartment for control of the water level in the reactor. As illustrated in Figure 1, wastewater was pumped into the first down-flow zone and then to the compartments one by one. The effluent from the last up-flow zone was discharged after the water level controlling tank. The biogas produced in each of the compartment was collected after the water lock. There were four equidistant sampling ports on one side of each up-flow zone for sampling solid and liquid. There was also an outlet at the bottom of each compartment for active sludge sampling. The reactor was wrapped with electrothermal wire and the temperature in the reactor was kept at 35 ± 1 °C by a temperature controller.

Figure 1

Schematic diagram of the anaerobic baffled reactor process: ① Feed tank; ② Peristaltic pump; ③ Baffle; ④ Gas meter; ⑤ Water lock; ⑥ Gas outlet; ⑦ Effluent; ⑧ Sampling port; ⑨ Solid outlet.

Figure 1

Schematic diagram of the anaerobic baffled reactor process: ① Feed tank; ② Peristaltic pump; ③ Baffle; ④ Gas meter; ⑤ Water lock; ⑥ Gas outlet; ⑦ Effluent; ⑧ Sampling port; ⑨ Solid outlet.

Close modal

The ABR was inoculated with 5.66 g mixed liquor volatile suspended solids (MLVSS) per liter of excess sludge that was collected from secondary clarifier of local sewage treatment plant (Harbin, China). Fed with synthetic wastewater through a peristaltic pump, the reactor was started up at a hydraulic retention time (HRT) of 2.0 d. An overloading of the reactor could lead to an inhibition of methanogens and consequently to failure of the start-up process. Therefore, normally the organic loading rate is low during the start-up period. The synthetic wastewater was brown sugar liquid with a chemical oxygen demand (COD) of about 2000 mg L−1. As illustrated in Table 1, the 101-day start-up period of the ABR was divided into five stages with an HRT of 2.0, 2.0, 1.0, 2.0 and 1.7 d, respectively. pH was adjusted to neutral through adding 3.36 g L−1 bicarbonate into influent during the last four stages.

Table 1

Operation stages of the ABR

StageTime (d)Influent COD (mg L−1)HRT (d)Organic load (kg COD m−3 d−1)NaHCO3 in influent (mg L−1)Influent pH
0–11 1,978 ± 200 2.0 1.00 ± 0.10 0 ± 0 6.60 ± 0.20 
12–30 2,145 ± 262 2.0 1.00 ± 0.13 3,360 ± 300 8.20 ± 0.30 
31–38 2,031 ± 180 1.0 2.00 ± 0.18 3,360 ± 200 8.20 ± 0.20 
39–75 4,246 ± 346 2.0 2.00 ± 0.17 3,360 ± 300 8.20 ± 0.30 
76–101 4,103 ± 254 1.7 2.40 ± 0.15 3,360 ± 200 8.20 ± 0.20 
StageTime (d)Influent COD (mg L−1)HRT (d)Organic load (kg COD m−3 d−1)NaHCO3 in influent (mg L−1)Influent pH
0–11 1,978 ± 200 2.0 1.00 ± 0.10 0 ± 0 6.60 ± 0.20 
12–30 2,145 ± 262 2.0 1.00 ± 0.13 3,360 ± 300 8.20 ± 0.30 
31–38 2,031 ± 180 1.0 2.00 ± 0.18 3,360 ± 200 8.20 ± 0.20 
39–75 4,246 ± 346 2.0 2.00 ± 0.17 3,360 ± 300 8.20 ± 0.30 
76–101 4,103 ± 254 1.7 2.40 ± 0.15 3,360 ± 200 8.20 ± 0.20 

Analytical methods

The influent and effluent of each compartment was sampled daily to test pH, COD and volatile fatty acid (VFA) contents. Analyses of COD and pH were carried out following standard methods (APHA 2005). The pH was determined using a pH meter (DELTA320, METTLER TOLEDO). Samples for VFA measurement were centrifuged at 13,000 rpm for 3 min using a Microfuge 16 Centrifuge (Beckman Coulter), and the supernatant was acidified by the addition of 25% H3PO4. VFA analysis was carried out using a gas chromatograph (SP6890, Shandong Lunan Instrument Factory, China) with injector, column and flame ionization detector (FID) temperatures at 210, 190 and 210 °C, respectively (Shi et al. 2019).

At the end of the start-up period, sludge samples were collected from the four compartments and the microbial communities were analyzed by high-throughput sequencing (Shanghai Bioengineering Technology Service Co., Ltd). The sequencing machine was an Illumina MiSeq PE (Illumina, USA). The 16S rRNA amplification primers for the Archaea were Arch349F (5′-GYGCASCAGKCGMGAAW-3′) and Arch806R (5′-GGACTA CVSGGGTATCTAAT-3′). Operational taxonomic units (OTUs) for each sequence were obtained from QIIME. Ribosomal Database Project (RDP) data were used to classify 16S rRNA gene sequences. The classification of functional microorganisms was based on studies by Müller et al. (2010), Stams et al. (2012) and Si et al. (2015).

Model implementation

The model, based on CA and ADM1, consists of the following main components (Figure 2).

  • (i)

    Compartment of ABR containing a liquid phase in contact with granular sludge. The bulk liquid compartment having a volume Vliq is continuously fed at a constant flow rate qin with a solution containing a number i of soluble substrates Sin,i. The bulk liquid is assumed to be completely mixed and continuously withdrawn from the reactor at the same flow rate (qin = qout). The bulk liquid volume is considerably larger than the granule volume, and therefore Vliq can be considered constant.

  • (ii)

    A field of concentrations of soluble compounds, based on a cartesian grid. The soluble compounds are the dissolved substrates and products involved in ADM1. This is the solution domain for all substrate components.

  • (iii)

    A set of biomass particles in different states. Each biomass particle (cell) is a separate differential state, represented by its position and type (species). Biomass particles are effectively independent of the substrate field. Each biomass particle is situated in one grid element with four possible directions in which it can spread.

Figure 2

The hypothetical granular sludge: Each biomass particle situated in one grid element had four directions to spread, () acidogens, () butyrate-degrading H2-producing acetogens, () propionate-degrading H2-producing acetogens, () hydrogenotrophic and () acetotrophic methanogens. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/wst.2021.113.

Figure 2

The hypothetical granular sludge: Each biomass particle situated in one grid element had four directions to spread, () acidogens, () butyrate-degrading H2-producing acetogens, () propionate-degrading H2-producing acetogens, () hydrogenotrophic and () acetotrophic methanogens. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/wst.2021.113.

Close modal

Domain and domain boundaries

The computational domain is only a small part of the whole system and includes both a small fraction of the bulk liquid and the granule (Figure 2(b)). The granule grows in a rectangular domain (box) with dimensions LX × LY (LX = 1,000 μm, LY = 1,000 μm). An ideal granule was assumed to have two-dimensional axial symmetry and therefore the computational domain was simplified to a quarter of the rectangular domain. The computational domain was divided into NX × NY regular grids (NX = 250, NY = 250) (Figure 2(c)). The mesh size of each grid was 2 μm × 2 μm, occupied by one cell and the solutes (Figure 2(d)).

Domain boundaries were treated as special conditions of biomass spreading. In the x-axis, biomass grew and spread in two directions:
formula
(1)
formula
(2)
where XL is left neighboring element, and XR is right neighboring element. When x is equal to 2, XL is defined as 1. When x is equal to NX − 1, XR is defined as NX.
Similarly, in the y-axis, biomass grew and spread in two directions:
formula
(3)
formula
(4)
where YD is the down neighboring element in the downward direction, and YU is the neighboring element in the upward direction. When y is equal to 2, YD is defined as 1. When y is equal to NY − 1, YU is defined as NY. Take the definition of boundary conditions in the y direction as an example. When the biomass was close to the domain boundary y = Ny (i.e., there is no neighboring element in the upward direction), then biomass growth and spreading stop.

Simulation of soluble components

Within each 2 × 2 μm grid element, solutes were transported in and out by diffusion (Fick's Law) to and from the neighboring elements. The diffusion distance was set at 2 μm. The time of biomass growth is usually much higher than that for diffusion of solutes (Batstone et al. 2006). Therefore, the mass balance involving diffusion can be assumed to be in a steady state each time (Batstone et al. 2006). The biochemical process was based on ADM1, expressed with 25 dynamic state variables and 20 biochemical rate equations as illustrated in Table S1 (Supplementary data). The mass balance in the ABR was described as follows (Stamatelatou et al. 2003):
formula
(5)
formula
(6)
formula
(7)
formula
(8)
where S1,i (kg COD m−3) is the concentration of soluble components in C1; qin and qout (m3 d−1) are the flow into and out of the compartment, respectively; Sin,i (kg COD m−3) is the input concentration of soluble components; Vliq (m3) is the liquid volume of the reactor; νi,j (kg COD m−3) is the stoichiometric coefficient; i is the serial number of the soluble components as illustrated in Table S1; and S2,i, S3,i and S4,i (kg COD m−3) are the concentration of the soluble components in C2, C3 and C4, respectively. Finally, ρj (kg COD m−3 d−1) is the specific kinetic rate for process j that represents the following 20 processes: disintegration; hydrolysis of carbohydrates, proteins and lipids; uptake of sugars, amino acids, long-chain fatty acid (LCFA), valerate, butyrate, propionate, acetate and hydrogen; and eight decays of biomass (Table S1).

Biomass growth and decay

Biomass involved in anaerobic digestion was assumed to be divided into the following six types: acidogens (Xsu), butyrate-degrading H2-producing acetogens (Xbu), propionate-degrading H2-producing acetogens (Xpro), hydrogenotrophic methanogens (Xh2), acetotrophic methanogens (Xac) and inert biomass (Xi). Xi resulted from the decay of the active biomass. The separation of active biomass into different categories was based on differences in metabolism. The growth and decay of each biomass type with time was described by an ordinary differential equation as follows:
formula
(9)
where Xi (kg COD m−3) is the mass concentration of Xsu, Xbu, Xpro, Xac and Xh2; rXi (kg COD m−3 d−1) is the biomass growth rate; and kdec,i (kg COD kg COD−1) is the decay rate of active biomass. Dynamic and stoichiometric parameters of biomass growth and decay were derived from Batstone et al. (2002) and Romli et al. (1995). An initial distribution of biomass particles at t = 0 must be defined. It was assumed that all initial biomass particles were randomly distributed in the core of granules. The initial coordinate of the granular core was LX = 500 μm and LY = 500 μm.

Biomass spreading

Spreading as a result of biomass growth occurs by occupying neighbor grid elements. In this two-dimensional model, biomass spreads to neighbor grids in four directions, up and down and left and right (Figure 2(d)). Biomass competence of neighbor grids was expressed as spreading probability (pi), determined by the ratio of each biomass type to the total biomass:
formula
(10)
where Xtotal (kg COD m−3) is the total biomass concentration.
An important interaction occurs between the syntrophic community simultaneously degrading higher organic acids to acetate and hydrogen or formate, and the methanogens that convert hydrogen and acetate to methane. Interspecific electron transfer was also considered when determining the spreading probability. From a reaction engineering and thermodynamic point of view, there is no difference between hydrogen and formate as electron carriers. Interspecific electron transfer was simplified to hydrogen mass transfer, which was described by Fick's Law as follows:
formula
(11)
where J (μmol μm−2 s−1) is interspecific hydrogen flux; D (m2 s−1) is the hydrogen diffusion coefficient, constant of 4.9 × 10−5 m2 s−1 with 298 K; A (μm−1) is the specific surface area of H2-producing acetogens; c1 (μmol m−2) is the hydrogen concentration on the surface of methanogens; c2 (μmol m−2) is the hydrogen concentration on the surface of H2-producing acetogens; and d (μm) is distance between H2-producing acetogens and methanogens, determined by biomass density. It was assumed that biomass density of methanogenic granular sludge is 1011–1012 cells ml−1 and rounded bacteria is 2 μm in diameter, consequently hydrogen diffusion distance is shorter than 0.5 μm (Stams et al. 2012). The hydrogen diffusion distance between species is less than the size of a grid element, therefore neighbor grid elements of hydrogenotrophic methanogens can be occupied by H2-producing acetogens. Batstone et al. (2006) showed that methanogen and H2-producing acetogen bacteria were observed as a physical co-culture within anaerobic granules and biofilms. Therefore the spreading probability of syntrophic bacteria can be calculated as follows:
formula
(12)
formula
(13)
formula
(14)
where pbu is the spreading probability of butyrate-degrading H2-producing acetogens; ppro is the spreading probability of propionate-degrading H2-producing acetogens; ph2 is the spreading probability of hydrogenotrophic methanogens; pbu,bu is the spreading probability of butyrate-degrading H2-producing acetogen reproduction; ppro,pro is the spreading probability of propionate-degrading H2-producing acetogen reproduction; ph2,h2 is the spreading probability of hydrogenotrophic methanogen reproduction; pbu,h2 is the syntrophic spreading probability between butyrate-degrading H2-producing acetogens and hydrogenotrophic methanogens; ppro,h2 is the syntrophic spreading probability between propionate-degrading H2-producing acetogens and hydrogenotrophic methanogens; ph2,bu is the probability that neighbor grid elements of butyrate-degrading H2-producing acetogens are occupied by hydrogenotrophic methanogens; and ph2,pro is the probability that neighbor grid elements of propionate-degrading H2-producing acetogens are occupied by hydrogenotrophic methanogens. Because the biomasses of Xaa, Xfa, Xc, Xpr, Xch, Xli, Xi and Xva are negligible in sugar fermentation, they are not considered in the spreading model.

Simulation

1-D model

In particular to describe VFA profiles in the ABR, the 1-D model was implemented in the software AQUASIM (Reichert 1998). The mixed reactor model provided by AQUASIM was used to simulate the mass transfer and conversion processes occurring in the bulk liquid of each compartment. The mixed reactors modeled for the compartments in series were connected by links and used to model the water and substance exchange among them. In this way, the biological conversion process in each compartment could be based on ADM1.

2-D model

The dynamics of soluble substrate concentrations relating to transport and reaction are rapid compared with the characteristic times of biomass growth, division, and spreading (Batstone et al. 2006). This observation allowed us to decouple the solution of biomass evolution from that of substrate diffusion reaction mass balances. Because of this time separation assumption, the algorithm for the modeling processes can be executed sequentially, as explained in the pseudocode below.

  • Set initial condition for biomass

  • Set initial condition for solutes in bulk liquid

  • do

  • 1-D model

    • 1.

      Biomass growth and decay (ADM1)

    • 2.

      Mass balances of solutes in the bulk liquid (ADM1)

    • 3.

      Physicochemical processes in the bulk liquid and gas phase (ADM1)

  • 2-D model

    • 4.

      Biomass spreading (CA)

    • 5.

      Modified biomass based on the ratio of each biomass type (CA + ADM1)

  • Time step, tt+ Δt

  • whilet < tend

First, random seed granules comprising 50 cells of each biomass type, with individual particle diameters of 2 μm were settled in compartments of the ABR. The bulk concentrations of monosaccharides are 2,000 mg COD L−1 at pH of 7. For biomass growth and decay, the biomass balances for each particle were then simulated based on ADM1. Solute concentrations were computed in the bulk liquid. At the same time step, physicochemical processes were also simulated in the bulk liquid. Many particles overlapped following the growth steps, and therefore it was necessary to simulate spreading according to the CA theory described above. After some biomass had reproduced and moved to the neighbor grid in the spreading step, the biomass dynamics in the time interval Δt were considered completed; that is, all new positions and distribution at t+ Δt were known. This time step lasted for 0.125 h. With this new biomass distribution, the proportion of each biomass type was estimated to modified biomass concentration in ADM1. We then computed new fields of solute concentrations by integrating CA to ADM1. Spreading ended (tend) once the biomass concentration had stabilized. According to CA rules, the inactive site would finally be occupied by active biomass. For simplification, biomass decay was not considered in the 2-D model.

Analysis of simulation outputs

To study the influence of interspecific competition on granule growth, the integrated model was used to simulate an ABR start-up. Two criteria were used to assess the simulations: (i) product concentration and (ii) biomass distribution. To provide a quantitative comparison, VFAs and methanogen biomass were both measured and predicted. It is generally agreed that VFAs accumulate because of biomass syntrophic growth and metabolism. At the end-stage of the start-up, the high throughput sequencing results of methanogens were used for model validation. Measured biomass of hydrogenotrophic and acetotrophic methanogens can provide direct evidence of interspecific competition. Comparison of predicted and measured biomass was also used as an assessment of simulation accuracy.

Interspecific competition in acidic and neutral conditions

Predictions of cell distribution within the granules in acidic and neutral conditions are shown in Figure 3. At the initial stage of reactor start-up, low pH of 4.5–4.6 substantially inhibited acetogens and methanogens. Acidogens grew in quantity and occupied a primary place in granules, whereas syntrophic communities were scattered, resulting in long distance hydrogen transfer and inhibition of syntrophic growth. The predicted ratios of acidogens of total biomass in C1, C2, C3 and C4 were 91.24%, 90.75%, 88.51% and 86.75%, respectively. In contrast, the ratios for hydrogenotrophic and acetotrophic methanogens in C1, C2, C3 and C4 were 8.22%, 8.67%, 10.64% and 12.59%, respectively (Figure 4). Acidogens were dominant in the reactor. When pH was adjusted to 6.7–7.3, the biomass of acetogens and methanogens increased. The central core of the granules was occupied by syntrophic communities, whereas acidogens were squeezed into the outer layer. A simulation in which acetogens and methanogens competed for the central region, resulted in a layered structure of granules. The predicted ratios of acidogens of the total biomass in C1, C2, C3 and C4 decreased to 40.32%, 37.12%, 30.45% and 29.73%, respectively, whereas the ratios of hydrogenotrophic and acetotrophic methanogens increased to 28.36%, 35.29%, 46.54% and 52.97%, respectively (Figure 4). The simulation results indicate that the biomass of acidogens decreased in the compartments under neutral pH conditions, whereas acetogens and methanogens were dominant in C4, which is a typical characteristic of bio-phase separation.

Figure 3

Spatial distribution of microflora in the ABR with acidic and neutral conditions: () Liquid phase, () Acidogens, () Butyrate-degrading H2-producing acetogens, () Propionate-degrading H2-producing acetogens, () Hydrogenotrophic methanogens, () Acetotrophic methanogens. (a) Compartment 1 (pH 4.5–4.6), (b) Compartment 2 (pH 4.5–4.6), (c) Compartment 3 (pH 4.5–4.6), (d) Compartment 4 (pH 4.5–4.6), (e) Compartment 1 (pH 6.7–7.3), (f) Compartment 2 (pH 6.7–7.3), (g) Compartment 3 (pH 6.7–7.3), (h) Compartment 4 (pH 6.7–7.3). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/wst.2021.113.

Figure 3

Spatial distribution of microflora in the ABR with acidic and neutral conditions: () Liquid phase, () Acidogens, () Butyrate-degrading H2-producing acetogens, () Propionate-degrading H2-producing acetogens, () Hydrogenotrophic methanogens, () Acetotrophic methanogens. (a) Compartment 1 (pH 4.5–4.6), (b) Compartment 2 (pH 4.5–4.6), (c) Compartment 3 (pH 4.5–4.6), (d) Compartment 4 (pH 4.5–4.6), (e) Compartment 1 (pH 6.7–7.3), (f) Compartment 2 (pH 6.7–7.3), (g) Compartment 3 (pH 6.7–7.3), (h) Compartment 4 (pH 6.7–7.3). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/wst.2021.113.

Close modal
Figure 4

Relative amount of microflora in the granule field simulated by the model: () Acidogens, () Butyrate-degrading H2-producing acetogens, () Propionate-degrading H2-producing acetogens, () Hydrogenotrophic methanogens, () Acetotrophic methanogens. (a) Compartment 1 (pH 4.5–4.6), (b) Compartment 2 (pH 4.5–4.6), (c) Compartment 3 (pH 4.5–4.6), (d) Compartment 4 (pH 4.5–4.6), (e) Compartment 1 (pH 6.7–7.3), (f) Compartment 2 (pH 6.7–7.3), (g) Compartment 3 (pH 6.7–7.3), (h) Compartment 4 (pH 6.7–7.3). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/wst.2021.113.

Figure 4

Relative amount of microflora in the granule field simulated by the model: () Acidogens, () Butyrate-degrading H2-producing acetogens, () Propionate-degrading H2-producing acetogens, () Hydrogenotrophic methanogens, () Acetotrophic methanogens. (a) Compartment 1 (pH 4.5–4.6), (b) Compartment 2 (pH 4.5–4.6), (c) Compartment 3 (pH 4.5–4.6), (d) Compartment 4 (pH 4.5–4.6), (e) Compartment 1 (pH 6.7–7.3), (f) Compartment 2 (pH 6.7–7.3), (g) Compartment 3 (pH 6.7–7.3), (h) Compartment 4 (pH 6.7–7.3). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/wst.2021.113.

Close modal

Relationship between interspecific competition and VFAs

The first 38-day start-up period of the ABR was divided into three stages as shown in Table 1, and the simulated VFAs of each stage are illustrated in Figure 5. The results show that acetate, propionate and butyrate were the primary liquid products in the ABR at Stage 1, with pH range 4.5 to 4.6. Predicted biomass of acetogens and methanogens was lower than 0.04 kg COD m−3. Following an increase in pH in the ABR from 6.7 to 7.3 in Stage 2, the concentration of propionate increased markedly, whereas the concentrations of butyrate and acetate decreased. The simulation showed that logarithmic growth in acetotrophic methanogens and butyrate-degrading H2-producing acetogens would occur from Day 11 to Day 18. These communities spread rapidly and competed for neighboring space, resulting in inhibition of propionate-degrading H2-producing acetogens. After Day 20, the biomass of butyrate-degrading H2-producing acetogens remained in a steady state, whereas the biomass of propionate-degrading H2-producing acetogens increased and propionate concentration subsequently declined sharply. The simulation results suggested that interspecific competition among acetogens had an effect on the butyrate and propionate concentrations.

Figure 5

Simulated outflow of VFAs and total biomass: () Butyrate measured, () Propionate measured, () Acetate measured, () Simulated results, () Predicted propionate-degrading H2-producing acetogens, () Predicted butyrate-degrading H2-producing acetogens, () Predicted acetotrophic methanogens. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/wst.2021.113.

Figure 5

Simulated outflow of VFAs and total biomass: () Butyrate measured, () Propionate measured, () Acetate measured, () Simulated results, () Predicted propionate-degrading H2-producing acetogens, () Predicted butyrate-degrading H2-producing acetogens, () Predicted acetotrophic methanogens. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/wst.2021.113.

Close modal

Interspecific competition between hydrogenotrophic and acetotrophic methanogens

High-throughput sequencing technology was used to analyze the distribution of methanogens in the four compartments of the ABR (Figure 6). The analysis results showed that nine families of Archaea (abundance larger than 1%): Fervidicoccaceae, Methanospirillaceae, Methanocellaceae, Methanomicrobiaceae, Methanoregulaceae, Methanococcaceae, Methanobacteriaceae, Methanosaetaceae and Methanosarcinaceae. Among these, Methanospirillaceae, Methanocellaceae, Methanomicrobiaceae, Methanoregulaceae, Methanococcaceae and Methanobacteriaceae are hydrogenotrophic methanogens, while Methanosaetaceae and Methanosarcinaceae are acetotrophic methanogens. The richness and diversity of the archaeal community in the ABR are shown in Table 2. Archaeal communities were similar in all four compartments. C2 had a relatively high richness, with the number of OTUs, ACE and Chao being 3,132, 8,534 and 6,207, respectively.

Table 2

Richness and diversity of archaeal community in the ABR

SampleSequence numbersRichness
Diversity index
Coverage (%)
OTUACEChaoShannonSimpson
Compartment 1 52,805 2,996 7,785 5,825 4.31 0.0514 96.96 
Compartment 2 57,357 3,132 8,534 6,207 4.32 0.0634 97.06 
Compartment 3 54,199 2,912 6,146 5,220 4.02 0.0898 97.39 
Compartment 4 56,850 2,962 7,465 5,646 4.04 0.1105 97.32 
SampleSequence numbersRichness
Diversity index
Coverage (%)
OTUACEChaoShannonSimpson
Compartment 1 52,805 2,996 7,785 5,825 4.31 0.0514 96.96 
Compartment 2 57,357 3,132 8,534 6,207 4.32 0.0634 97.06 
Compartment 3 54,199 2,912 6,146 5,220 4.02 0.0898 97.39 
Compartment 4 56,850 2,962 7,465 5,646 4.04 0.1105 97.32 
Figure 6

Relative amount of Archaea at the family level, and ‘others’ includes the Archaea that were less than 1% relative amount of abundance. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/wst.2021.113.

Figure 6

Relative amount of Archaea at the family level, and ‘others’ includes the Archaea that were less than 1% relative amount of abundance. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/wst.2021.113.

Close modal

The ratios of hydrogenotrophic methanogens in C1, C2, C3 and C4 were 57.43%, 36.09%, 35.16% and 34.02%, respectively. Acetotrophic methanogens were 42.57%, 63.91, 64.84 and 65.98%, respectively (Figure 7). Simulation results showed a decline in the hydrogenotrophic methanogens through the compartments. In C1, hydrogenotrophic methanogens were dominant. This suggests that hydrogenotrophic methanogens were more competitive than acetotrophic methanogens with an adequate supply of substrate. It has been reported that Methanomicrobiaceae, Methanobacteriaceae, Methanosaetaceae and Methanosarcinaceae were the most common methanogens in the ABR (Ziganshin et al. 2019). Most of the methanogens, particularly acetotrophic methanogens, were predominantly distributed in the last compartments, where VFA concentration considerably decreased and the highest biogas production was observed (Tran et al. 2017). Conversely, growth of hydrogenotrophic methanogens was inhibited by low hydrogen partial pressure, resulting in acetotrophic methanogens being dominant in the last compartment. The results indicated that interspecific competition affected the community structure of methanogens.

Figure 7

Comparison between sequencing results and predicted data of methanogens biomass in ABR: (a) ratio of hydrogenotrophic methanogens and (b) ratio of acetotrophic methanogens.

Figure 7

Comparison between sequencing results and predicted data of methanogens biomass in ABR: (a) ratio of hydrogenotrophic methanogens and (b) ratio of acetotrophic methanogens.

Close modal

The layered structure of the microbial community

The results of the simulations indicated that, contrary to the syntrophic model, interspecific competition is an important mechanism for the formation of layered microbial community structures. In addition, syntrophic communities compete for core space in favor of interspecific electron transfer. The layered structure of the microbial community supports the findings of Arching et al. (1993) and Lens et al. (1995). It was assumed that the inner layer mainly consists of methanogens that act as the nucleation centers necessary for the initiation of granule development. In the current simulations, the inner layer structure also resulted from interspecific competition among acidogens, acetogens and methanogens. Interspecific competition generally occurs among acetogens and methanogens, and has an important effect on the syntrophic community. However, the effect of interspecific competition is reversible. The results of this study show that the cellular automaton model can simulate a dynamic development of granules under varying environmental conditions and it also produced a large variety of distinct morphologies in response to changes in growth conditions.

In this paper we describe the development of an integrated model combining CA and ADM1 for dynamic modeling of interspecific competition in the microbial community of an ABR. Predictions of VFAs and methanogenic biomass were quantitatively and critically compared with experimental data obtained from an ABR start-up. Syntrophic co-location of hydrogen producers and consumers and the multi-layer structure of granules were clearly observed in the model. Butyrate-degrading acetogens grew in quantity, resulting in the inhibition of growth of propionate-degrading acetogens. In syntrophic communities, hydrogenotrophic methanogens showed a stronger competitive advantage than acetotrophic methanogens. The biomass distribution of methanogens was predicted by the model.

Supplementary data for this work can be found in online version of the paper.

The authors gratefully acknowledge the National Science and Technology Major Project, China (grant number 2018ZX07601001), Fundamental Research Funds for Colleges and Universities in Liaoning Province, China (grant number lnjc202012) and Fundamental Research Funds for Colleges and Universities in Liaoning Province, China (grant number lnqn202020).

All relevant data are included in the paper or its Supplementary Information.

APHA, AWWA, WPCF
2005
Standard Methods for the Examination of Water and Wastewater
, 21st edn.
Water Research, Washington, DC
.
Arching
B. K.
Schmidt
J. E.
Winther-Nielsen
M.
Macario
A. J. L.
de Macario
E. C.
1993
Effect of medium composition and sludge removal on the production, composition and architecture of thermophilic (551c) acetate-utilizing granules from an upflow anaerobic sludge blanket reactor
.
Applied and Environmental Microbiology
59
,
2538
2545
.
Baeten
J. E.
Batstone
D. J.
Schraa
O. J.
van Loosdrecht
M. C. M.
Volcke
E. I. P.
2019
Modelling anaerobic, aerobic and partial nitritation-anammox granular sludge reactors – A review
.
Water Research
149
,
322
341
.
Batstone
D. J.
Keller
J.
Angelidaki
I.
Kalyuzhnyi
S. V.
Pavlostathis
S. G.
Rozzi
A.
Sanders
W. T.
Siegrist
H.
Vavilin
V. A.
2002
Anaerobic Digestion Model No. 1 (ADM1)
.
IWA Publishing
,
London
.
Batstone
D. J.
Keller
J.
Blackall
L. L.
2004
The influence of substrate kinetics on the microbial community structure in granular anaerobic biomass
.
Water Research
38
,
1390
1404
.
Batstone
D. J.
Picioreanu
C.
van Loosdrecht
M. C. M.
2006
Multidimensional modelling to investigate interspecies hydrogen transfer in anaerobic biofilms
.
Water Research
40
,
3099
3108
.
Batstone
D. J.
Puyol
D.
Flores-Alsina
X.
Rodríguez
J.
2015
Mathematical modelling of anaerobic digestion processes: applications and future needs
.
Reviews in Environmental Science and Bio/Technology
14
(
4
),
595
613
.
Cook
S. M.
Skerlos
S. J.
Raskin
L.
Love
N. G.
2017
A stability assessment tool for anaerobic codigestion
.
Water Research
112
,
19
28
.
Fang
H. H. P.
2000
Microbial distribution in UASB granules and its resulting effects
.
Water Science and Technology
42
,
201
208
.
Fang
H. H. P.
Ho-Kwong
C.
Yu-You
L.
1995
Effect of degradation kinetics on the microstructure of anaerobic biogranules
.
Water Science and Technology
32
(
8
),
165
172
.
Feldman
H.
Flores-Alsina
X.
Ramin
P.
Kjellberg
K.
Jeppsson
U.
Batstone
D. J.
Gernaey
K. V.
2017
Modelling an industrial anaerobic granular reactor using a multi-scale approach
.
Water Research
126
,
488
500
.
Gagliano
M. C.
Ismail
S. B.
Stams
A. J. M.
Plugge
C. M.
Temmink
H.
Van Lier
J. B.
2017
Biofilm formation and granule properties in anaerobic digestion at high salinity
.
Water Research
121
,
61
71
.
Grotenhuis
J. T. C.
Lier
J. B. V.
Plugge
C. M.
Stams
A. J. M.
Zehnder
A. J. B.
1991
Effect of ethylene glycol-bis(β-aminoethyl ether)-N,N-tetraacetic acid (EGTA) on stability and activity of methanogenic granular sludge
.
Applied Microbiology and Biotechnology
36
(
1
),
109
114
.
Guiot
S. R.
Pauss
A.
Costerton
J. W.
1992
A structured model of the anaerobic granules consortium
.
Water Science and Technology
25
,
1
10
.
Hulshoff Pol
L. W.
de Castro Lopes
S. I.
Lettinga
G.
Lens
P. N. L.
2004
Anaerobic sludge granulation
.
Water Research
38
,
1376
1389
.
Jeppsson
U.
Alex
J.
Batstone
D. J.
Benedetti
L.
Comas
J.
Copp
J. B.
Corominas
L.
Flores-Alsina
X.
Gernaey
K. V.
Nopens
I.
Pons
M.-N.
Rodriguez-Roda
I.
Rosen
C.
Steyer
J.-P.
Vanrolleghem
P. A.
Volcke
E. I. P.
Vrecko
D.
2016
Benchmark simulation models, quo vadis?
Water Science & Technology
68
(
1
),
1
15
.
Laspidou
C. S.
Kungolos
A.
Samaras
P.
2010
Cellular-automata and individual-based approaches for the modeling of biofilm structures: pros and cons
.
Desalination
250
(
1
),
390
394
.
Lens
P.
de Beer
D.
Cronenberg
C.
Ottengraf
S.
Verstraete
W.
1995
The use of microsensors to determine distributions in UASB aggregates
.
Water Science and Technology
31
,
273
280
.
MacLeod
F. A.
Guiot
S. R.
Costerton
J. W.
1990
Layered structure of bacterial aggregates produced in an upflow anaerobic sludge bed and filter reactor
.
Applied and Environmental Microbiology
56
,
1598
1607
.
Manukyan
L.
Montandon
S. A.
Fofonjka
A.
Smirnov
S.
Milinkovitch
M. C.
2017
A living mesoscopic cellular automaton made of skin scales
.
Nature
544
(
7649
),
173
179
.
Müller
N.
Worm
P.
Schink
B.
Stams
A. J. M.
Plugge
C. M.
2010
Syntrophic butyrate and propionate oxidation processes: from genomes to reaction mechanisms
.
Environmental Microbiology Reports
2
(
4
),
489
499
.
Odriozola
M.
López
I.
Borzacconi
L.
2016
Modeling granule development and reactor performance on anaerobic granular sludge reactors
.
Journal of Environmental Chemical Engineering
4
,
1615
1628
.
Pauss
A.
Samson
R.
Guiot
S.
1990
Thermodynamic evidence of trophic microniches in methanogenic granular sludge-bed reactors
.
Applied Microbiology and Biotechnology
33
,
88
92
.
Picioreanu
C.
Loosdrecht
M. C. M. V.
Heijnen
J. J.
1998
A new combined differential-discrete cellular automaton approach for biofilm modeling: application for growth in gel beads
.
Biotechnology and Bioengineering
57
(
6
),
718
731
.
Picioreanu
C.
Batstone
D. J.
van Loosdrecht
M. C. M.
2005
Multidimensional modelling of anaerobic granules
.
Water Science and Technology
52
,
501
507
.
Reichert
P.
1998
AQUASIM 2.0-Computer Program for the Identification and Simulation of Aquatic Systems
.
Version 2.0
.
EAWAG
,
Dubendorf
,
Switzerland
.
Romli
M.
Keller
J.
Lee
P. J.
Greenfield
P. F.
1995
Model prediction and verification of a two-stage high-rate anaerobic wastewater treatment system subjected to shock loads
.
Process Safety and Environmental Protection
73
(
2
),
151
154
.
Sanchez
J. M.
Arijo
S.
Munoz
M. A.
Morinigo
M. A.
Borrego
J. J.
1994
Microbial colonization of different support materials used to enhance the methanogenic process
.
Applied Microbiology and Biotechnology
41
(
4
),
480
486
.
Schmidt
J. E.
Ahring
B. K.
1996
Granular sludge formation in upflow anaerobic sludge blanket (UASB) reactors
.
Biotechnology and Bioengineering
49
(
3
),
229
246
.
Si
B.
Liu
Z.
Zhang
Y.
Li
J.
Xing
X. H.
Li
B.
Duan
N.
Lu
H.
2015
Effect of reaction mode on biohydrogen production and its microbial diversity
.
International Journal of Hydrogen Energy
40
(
8
),
3191
3200
.
Stamatelatou
K.
Lokshina
L.
Vavilin
V.
Lyberatos
G.
2003
Performance of a glucose fed periodic anaerobic baffled reactor under increasing organic loading conditions: 2. Model prediction
.
Bioresource Technology
88
,
137
142
.
Stams
A. J. M.
Worm
P.
Sousa
D. Z.
Alves
M. M.
Plugge
C. M.
2012
Syntrophic degradation of fatty acids by methanogenic communities
.
Microbial Technologies in Advanced Biofuels Production, Chapter
8
,
127
142
.
Tang
Y.
Valocchi
A. J.
2013
An improved cellular automaton method to model multispecies biofilms
.
Water Research
47
(
15
),
5729
5742
.
Tran
P. T.
Watari
T.
Hirakata
Y.
Nguyen
T.
Hatamoto
M.
Tanikawa
D.
Syutsubo
K.
Nguyen
T.
Fukuda
M.
Nguyen
L.
Yamaguchi
T.
2017
Anaerobic baffled reactor in treatment of natural rubber processing wastewater: reactor performance and analysis of microbial community
.
Journal of Water and Environment Technology
15
,
241
251
.
Wu
W.
Jain
M. K.
Zeikus
J. G.
1996
Formation of fatty acid-degrading, anaerobic granules by defined species
.
Applied and Environmental Microbiology
62
(
6
),
2037
2044
.
Young
J. C.
McCarty
P. L.
1969
The anaerobic filter for waste treatment
.
Journal of the Water Pollution Control Federation
41R
,
160
173
.
Ziganshin
A. M.
Wintsche
B.
Seifert
J.
Carstensen
M.
Born
J.
Kleinsteuber
S.
2019
Spatial separation of metabolic stages in a tube anaerobic baffled reactor: reactor performance and microbial community dynamics
.
Applied Microbiology and Biotechnology
103
(
9
),
3915
3929
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC-ND 4.0), which permits copying and redistribution for non-commercial purposes with no derivatives, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Supplementary data