Modeling interspecific competition of the microbial community during anaerobic digestion based on cellular automata and ADM1

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.

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.

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⁻¹. 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⁻¹ bicarbonate into influent during the last four stages.

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% H₃PO₄. 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).

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 V_liq is continuously fed at a constant flow rate q_in with a solution containing a number i of soluble substrates S_in,i. The bulk liquid is assumed to be completely mixed and continuously withdrawn from the reactor at the same flow rate (q_in = q_out). The bulk liquid volume is considerably larger than the granule volume, and therefore V_liq 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.

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 L_X × L_Y (L_X = 1,000 μm, L_Y = 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 N_X × N_Y regular grids (N_X = 250, N_Y = 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)).

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.

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, t → t+ Δt

• whilet < t_end

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⁻¹ 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 (t_end) 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.

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.

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.

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⁻³. 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 H₂-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 H₂-producing acetogens. After Day 20, the biomass of butyrate-degrading H₂-producing acetogens remained in a steady state, whereas the biomass of propionate-degrading H₂-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.

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.

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.

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.